# Generation de nombres aleatoires autocorreles

Marc Badel

► To cite this version:

Marc Badel. Generation de nombres aleatoires autocorreles. [Rapport de recherche] IRIA-RR-265, IRIA. 1977, pp.25. hal-04716471

HAL Id: hal-04716471

https://inria.hal.science/hal-04716471v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dw-03641 RR265

RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

GÉNÉRATION DE
NOMBRES ALÉATOIRES
AUTOCORRÉLÉS

Marc BADEL

Rapport de Recherche N° 265

Novembre 1977

esp.
GÉNÉRATION DE NOMBRES ALÉATOIRES AUTOCORRELES

Marc Badel

IRIA/LABORIA

Résumé :

De nombreuses simulations sur ordinateur nécessitent la génération d'échantillons artificiels destinés à représenter des phénomènes aléatoires à l'entrée d'un système physique. Le plus souvent, ces échantillons artificiels sont indépendants et de nombreux algorithmes de génération ont déjà été proposés pour ces cas. Il devient cependant parfois nécessaire de s'affranchir de cette indépendance. Un cas particulier important de non indépendance d'un échantillon est l'autocorrélation des grandeurs successives obtenues. Nous traitons ici le problème de la génération d'échantillons de distribution et d'autocorrélation donnés. Nous exposons d'abord un algorithme existant valable dans un cas particulier. Puis nous étendons cet algorithme à d'autres formes de distribution plus générales.

Abstract :

Many computer simulations need to generate artificial samples for representing random phenome-
na. These samples are usually assumed to be identically and independtly distributed. Numerous
algorithms have been presented which solve the problem in this specific case. However it is some-
times absolutely necessary to relax the independence assumption, for example if the sample has
to be autocorrelated. This paper is dealing with this problem of generating an artificial sample the
distribution and autocorrelation functions of which are given. We first present an existing algori-
thm available only in a particular case. Then we extend this algorithm to more general distribu-
tions.
# I - INTRODUCTION

La simulation sur ordinateur est l'une des techniques utilisées pour reproduire le fonctionnement d'un système modélisé : études de gestions de stocks, de réseaux de files d'attente, de chaînes de production ; ou l'évolution dans le temps d'un phénomène physique : études de désintégrations de corps radio-actifs, de populations. On pourrait multiplier les exemples.

Etant donnée l'approche probabiliste qui est souvent prise pour l'étude de ces problèmes, on se trouve confronté au paradoxe d'avoir à engendrer du hasard (des échantillons de variables aléatoires) à partir de mécanismes déterministes (un ordinateur).

Notons au passage que des méthodes, dites de Monte Carlo, permettent de résoudre des problèmes où, a priori, le hasard n'intervient pas. C'est le cas, par exemple, de l'intégration de surfaces ou de volumes compliqués par des méthodes d'échantillonnage statistique.

Cette nécessité de reproduire artificiellement le hasard a donc posé très tôt le problème de la génération de suites de nombres aléatoires. Des algorithmes ont été proposés pour générer de telles suites. On pourra trouver des précisions à ce sujet dans [5, 7].

Deux qualités de base sont demandées aux suites de nombres au hasard (nous parlerons désormais de suites pseûdo-aléatoires pour indiquer le caractère déterministe des algorithmes de génération) :

1. Distribution identique des nombres générés,

2. Indépendance stochastique de ces nombres.
Suivant les problèmes, la première propriété peut seule suffire (intégrations de surfaces ou de volumes), si aucune hypothèse d'indépendance n'est faite.

En ce qui concerne la deuxième propriété, la plupart des travaux ont recherché les manières de générer des suites de nombres indépendants et de s'assurer par des tests de cette indépendance [4, 6].

Cette propriété d'indépendance ou de non-indépendance est primordiale lorsqu'on s'intéresse à des phénomènes aléatoires dépendants du temps (processus stochastiques). Le cas de l'indépendance étant plus ou moins résolu (ce que nous supposerons), nous nous intéressons, dans ce papier, à la génération de suites de nombres pseudo-aléatoires arbitrairement autocorrélés.

Les deux propriétés que nous demandons aux suites que nous étudions sont donc les suivantes :

1 - nombres identiquement distribués suivant une loi de distribution donnée,

2 - nombres autocorrélés suivant une fonction d'autocorrélation donnée.

Une méthode pour engendrer de telles suites consiste à se ramener à des suites indépendantes. Nous nous limitons au cas d'autocorrélations à décroissance géométrique. Nous rappelons d'abord un algorithme existant, valable dans un cas particulier : génération d'échantillons autocorrélés de distribution exponentielle négative. Puis nous généralisons cet algorithme à d'autres formes de distributions usuelles. Ces algorithmes nous conduisent dans une phase finale à étendre la définition de certaines distributions connues et à proposer de nouvelles méthodes de génération d'échantillons indépendants pour ces cas.
## II - PRESENTATION DU PROBLEME

### II.1 - Définition d'un échantillon autocorrélé

Nous nous proposons de générer un échantillon $x_0, x_1, \dots, x_n, \dots$ d'une suite de variables aléatoires $X_0, X_1, \dots, X_n, \dots$ autocorrélées identiquement distribuées, c'est-à-dire :

- $X_0, X_1, \dots, X_n, \dots$ ont toutes la même distribution ou bien, plus simplement, les mêmes moyennes et variances si nous ne nous intéressons qu'aux premiers moments.

- La corrélation entre $X_n$ et $X_{n+1}$ a une valeur donnée $\rho$ indépendante de l'ordre n de la suite, c'est-à-dire :

$$ \text{Corr}(X_n, X_{n+1}) = \rho \text{ ou } \text{Cov}(X_n, X_{n+1}) = \rho V(X_n) $$

### II.2 - Génération d'un échantillon autocorrélé distribué exponentiellement

Un algorithme dû à Jacobs et Lewis [3] consiste à tirer $x_{n+1}$ à partir de $x_n$ avec

$$ (1) \quad x_{n+1} = \begin{cases} \rho x_n & \text{avec une probabilité } \rho \\ \rho x_n + u_n & \text{avec une probabilité } (1-\rho) \end{cases} $$

où $u_n$ est une suite de réalisations indépendantes de la variable aléatoire U distribuée exponentiellement avec :

$$ E(U) = E(X_n) $$

Nous avons :

$$ \text{Corr}(X_n, X_{n+1}) = \rho \quad \text{et} \quad \text{Corr}(X_n, X_{n+p}) = \rho^p $$

La démonstration de ces propriétés peut résulter du cas général développé ci-après.
## II.3 - Vers une généralisation de l'algorithme au cas non-exponentiel

L'algorithme (1) permet ainsi, dans le cas de suites autocorrélées de variables aléatoires distribuées exponentiellement, de se ramener à la génération d'une suite de variables aléatoires indépendantes et identiquement distribuées. Nous nous proposons de généraliser cet algorithme de façon à générer :

- Une suite $X_n$ de variables aléatoires autocorrélées, d'espérance $\xi$ et de variance $\gamma$ données

ou

- Une suite $X_n$ de variables aléatoires autocorrélées et identiquement distribuées de densité $f(x)$ donnée.

L'algorithme général que nous nous proposons de mettre en oeuvre à la forme :

$$ (2) \qquad x_{n+1} = \rho x_n + y_n $$

où $y_n$ est une suite de réalisations indépendantes et identiquement distribuées d'une variable aléatoire Y dont la fonction de densité $g(x)$ doit être déterminée de façon à satisfaire aux conditions imposées. Nous allons traiter successivement le problème des suites autocorrélées d'espérance et de variance données, puis celui des suites autocorrélées de distribution donnée.

## III - GENERATION D'UN ECHANTILLON AUTOCORRELE D'ESPERANCE ET DE VARIANCE DONNEES

### III.1 - Généralités

Dans de nombreux problèmes pratiques, nous nous intéressons peu à la forme des distributions mises en jeu, mais seulement à leurs premiers moments, habituellement leur moyenne et leur variance. Nous nous proposons d'utiliser l'algorithme (2) pour générer un échantillon $\{x_n\}$ d'une suite de variables aléatoires $X_0, X_1, \dots, X_n$ d'espérance $\xi$ et de variance $\gamma$ données (la méthode mise en oeuvre se généralise immédiatement à un nombre de moments plus élevé). Nous devons imposer pour cela des conditions convenables à $\rho$ et aux premiers moments de la variable aléatoire Y.
### III.2 - Conditions à imposer à la distribution des $y_n$

Puisque nous ne nous intéressons qu'aux deux premiers moments des $X_n$, les conditions à imposer portent également sur les deux premiers moments de Y. La relation (2) donne :

$$E(X_{n+1}) = \rho E(X_n) + E(Y)$$

et comme nous imposons $E(X_{n+1}) = E(X_n) = \xi$, nous avons :

$$(3) \quad E(Y) = (1-\rho) E(X_n) = (1-\rho) \xi$$

D'autre part, comme $X_n$ et $Y$ sont des variables indépendantes :

$$V(X_{n+1}) = \rho^2 V(X_n) + V(Y)$$

et comme par hypothèse $V(X_n+1) = V(X_n) = \gamma$, $V(Y)$ est donné par :

$$(4) \quad V(Y) = (1-\rho^2) V(X_n) = (1-\rho^2) \gamma$$

Y peut être une variable aléatoire de distribution quelconque, pourvu que ses deux premiers moments satisfassent (3) et (4).

### III.3 - Calcul de l'autocorrélation de la suite obtenue

L'équation (2) entraîne :

$$\text{Cov}(X_{n+1}, X_n) = \rho \text{Cov}(X_n, X_n) + \text{Cov}(Y, X_n)$$

Or $X_n$ et $Y$ sont des variables indépendantes, d'où

$$\text{Cov}(Y, X_n) = 0$$
D'autre part, $\text{Cov}(X_n, X_n) = V(X_n)$ d'où

$$
(5) \quad \text{Cov}(X_{n+1}, X_n) = \rho \text{V}(X_n) \text{ ou bien}
$$

$$
(6) \quad \operatorname{Corr}(X_{n+1}, X_n) = \rho
$$

Le coefficient ρ doit donc être choisi égal à l'autocorrélation
désirée entre deux termes consécutifs de la suite à générer.

Calculons l'autocorrélation d'ordre p soit $\text{Corr}(X_{n+p}, X_n)$. L'équation (2) à l'ordre p donne :

$$
x_{n+p+1} = \rho x_{n+p} + y_{n+p} \text{ d'où}
$$

$$
\text{Cov}(X_{n+p+1}, X_n) = \rho \text{Cov}(X_{n+p}, X_n) + \text{Cov}(Y, X_n). \text{ Comme } Y \text{ et } X_n \text{ sont des variables indépendantes, } \text{Cov}(Y, X_n) = 0 \text{ d'où :}
$$

$$
\mathrm{Cov}(X_{n+p+1}, X_n) = \rho \mathrm{Cov}(X_{n+p}, X_n) \text{ d'où l'on tire :}
$$

$$
(7) \quad \mathrm{Cov}(X_{n+p}, X_n) = \rho^p \mathrm{V}(X_n) \text{ ou bien}
$$

$$
(8) \quad \text{Corr}(X_{n+p}, X_n) = \rho^p
$$

L'autocorrélation de la suite des $X_n$ générée est donc à décroissance géométrique et $\rho$ est le coefficient d'autocorrélation à l'ordre 1.

### III.4 - Convergence de la distribution des $X_n$ lorsque n tend vers l'infini

Notons $g(x)$ la fonction de densité de $Y$, $g*(s)$ sa transformée de Laplace, $f_n(x)$ la fonction de densité de $X_n$, $f_n*(s)$ sa transformée de Laplace.
Ecrivons d'abord la relation entre les distributions de $X_{n+1}$, $X_n$ et Y. Ces variables étant liées en raison de (2), la fonction de densité de $\rho_{X_n}$ est $\frac{1}{\rho} f_n(\frac{X_n}{\rho})$ et sa transformée de Laplace est $f_n^*(\rho s)$. D'après (2), la distribution de $X_{n+1}$ est le produit de convolution des distributions de $\rho_{X_n}$ et de Y, puisque ces variables sont indépendantes. En passant aux transformées de Laplace, il vient :

$$ (9) \quad f_{n+1}^*(s) = f_n^*(\rho s) \cdot g^*(s) \quad \text{ce qui entraîne :} $$

$$ f_n^*(s) = f_0^*(\rho^n s) \prod_{i=1}^{n} g^*(\rho^{i-1} s) $$

Etudions $f_n^*(s)$ quand $n$ tend vers l'infini. Comme nous avons supposé $0 < \rho < 1$, $\rho^n s$ tend vers 0 et $f_0^*(\rho^n s)$ tend vers $f_0^*(0) = 1$. Montrons que le produit infini

$$ P(s) = \prod_{i=1}^{\infty} g^*(\rho^{i-1}s) \text{ existe} $$

Nous avons $\text{Log}P(s) = \sum_{n=0}^{\infty} \text{Log } g^*(\rho^{i-1}s)$. Or, au voisinage de 0, $g^*(u) = g^*(0) + u \frac{dg^*(0)}{ds}(0) + O(u^2)$ avec $g^*(0) = 1$ et $\frac{dg^*(0)}{ds} = -E(Y)$ d'où si $i$ tend vers l'infini :

$\text{Log } g^*(\rho^{i-1}s) \sim -E(Y)\rho^{i-1}s$, terme général d'une série géométrique convergente.

La série donnant $\text{Log } P(s)$ converge, donc le produit infini $P(s)$ existe et nous pouvons écrire :

$$ (10) \quad f_{\infty}^*(s) = \prod_{i=1}^{\infty} g^*(\rho^{i-1}s) $$

La distribution des $X_n$ tend donc vers une distribution limite lorsque $n$ tend vers l'infini. Cette distribution ne dépend pas de la distribution initiale suivie par $X_0$, mais seulement de la distribution de Y.
D'autre part, d'après (10), cette distribution vérifie :

$$f_{\infty}^{*}(s) = f_{\infty}^{*}(\rho s) g^{*}(s)$$

ce qui montre que l'égalité (9) reste vraie à la limite. Prendre $f_0(x) = f_{\infty}(x)$ comme distribution initiale de $X_0$ conduit à la génération de suites autocorrélées et identiquement distribuées. Nous allons aborder maintenant ce problème.

# IV - GENERATION DE SUITES AUTOCORRELEES IDENTIQUEMENT DISTRIBUEES ET DE DISTRIBUTION DONNEE

## IV.1 - Présentation du problème

Le problème est de déterminer une distribution g pour Y, telle qu'en appliquant l'algorithme (2), tous les $X_n$ de la suite générée suivent une distribution donnée. Ce cas est un cas particulier du cas précédent (génération d'une suite d'espérance et de variance données) et toutes les formules données précédemment s'appliquent. En particulier, les distributions de $X_n$ et $X_{n+1}$ sont liées entre elles par la formule (9) et si f est la fonction de densité imposée aux $X_n$, la relation (9) donne :

$$(11) \quad f^{*}(s) = f^{*}(\rho s) \quad \hat{g}^{*}(s) \text{ ou}$$

$$(12) \quad g^{*}(s) = \frac{f^{*}(s)}{f^{*}(\rho s)}$$

équation donnant la transformée de Laplace de la fonction de densité de Y dès que la fonction de densité f des $X_n$ est imposée.

## IV.2 - Conditions pour que le problème soit soluble

La génération par (2) d'un échantillon $\{x_n\}$ d'une suite de variables $X_n$ de densité f est possible si et seulement si la formule (12) aboutit, après la prise de la transformée de Laplace inverse de g, à une fonction g qui soit une fonction de densité. Une condition nécessaire et suffisante peut être
trouvée dans [2, page 588] dont nous rappellerons les principaux points (les résultats sur les fonctions caractéristiques donnés dans [2] se généralisant aux transformées de Laplace).

$\frac{f^*(s)}{f^*(\rho s)}$ (avec $0 < \rho < 1$) est la transformée de Laplace d'une fonction de densité g si et seulement si la distribution définie par f appartient à une classe de distributions particulières : la classe L.

Par définition, une distribution appartient à la classe L si et seulement si elle est la distribution limite d'une variable aléatoire :

$S_n^* = (S_n - b_n)/a_n$, $a_n$ et $b_n$ étant des suites de coefficients vérifiant

$$ \lim_{n \to \infty} a_n = \infty \quad \text{et} \quad \lim_{n \to \infty} \frac{a_n + 1}{a_n} = 1 \quad \text{et où} $$

$$ S_n = U_1 + \dots + U_n $$

$U_1, \dots, U_n$ étant une suite de variables aléatoires mutuellement indépendantes.

Feller donne dans [2] un ensemble de conditions nécessaires et suffisantes d'appartenance d'une distribution à la classe L. La condition la plus importante est l'infinie divisibilité de la distribution donnée.

## IV.3 - Application à quelques distributions particulières

### IV.3.1 - Distribution gaussienne

Le cas est résolu dès le paragraphe III. Il suffit de remarquer dans (2) que si $X_n$ et $Y$ suivent des distributions gaussiennes, $X_{n+1} = \rho X_n + Y$ suit également une distribution gaussienne. L'algorithme (2) s'applique donc dans ce cas en tirant les $y_n$ suivant une distribution gaussienne dont la moyenne et la variance sont données par (3) et (4).
### IV.3.2 - Distribution exponentielle négative

Nous retrouvons l'algorithme (1) dans ce cas. Soit $1/\lambda$ l'espérance des $X_n$, soit :

$$f(x) = \lambda e^{-\lambda x} \text{. Nous avons :}$$

$$f^*(s) = \frac{\lambda}{\lambda+s} \text{ et en appliquant (12)}$$

$$g^*(s) = \frac{\frac{\lambda}{\lambda+s}}{\frac{\lambda}{\lambda+s} + \rho s} = \frac{\lambda + \rho s}{\lambda + s} = \rho + (1-\rho) \frac{\lambda}{\lambda+s}$$

$g^*(s)$ est la transformée de Laplace de la densité d'une variable Y telle que :

- $Y=0$ avec une probabilité $\rho$

- Avec une probabilité de $(1-\rho)$, Y suit une loi exponentielle négative d'espérance $1/\lambda$.

Nous retrouvons l'algorithme (1).

### IV.3.3 - Distribution d'Erlang d'ordre p

Notons $\frac{p}{\lambda}$ la moyenne de la distribution. Une telle distribution peut être interprétée comme la distribution de la somme de p variables aléatoires indépendantes de distributions exponentielles négatives et d'espérances $\frac{1}{\lambda}$, soit :

$$f^*(s) = \left(\frac{\lambda}{\lambda+s}\right)^p \text{. D'où en appliquant (12)}$$

$$g^*(s) = \left(\frac{\lambda + \rho s}{\lambda + s}\right)^p = [\rho + (1-\rho) \frac{\lambda}{\lambda+s}]^p$$

La distribution g de Y est la distribution de la somme de p variables aléatoires indépendantes égales à 0 avec une probabilité $\rho$ et tirées avec une probabilité de $(1-\rho)$ suivant une loi exponentielle négative d'espérance $1/\lambda$.
D'où un algorithme de tirage de $y_n$ pour appliquer (2) dans ce cas.

### IV.3.4 - Distributions de Cox [1]

Par définition, les distributions de Cox sont les distributions dont la densité est à transformée de Laplace rationnelle (les distributions exponentielles et d'Erlang en sont des cas particuliers). D'après (25), si $f^*(s)$ est une fonction rationnelle de s, $g^*(s)$ est également rationnelle en s. Si (12) définit effectivement une distribution de probabilité, c'est-à-dire si la condition donnée en IV.2 (notamment l'infinie divisibilité de f) est vérifiée, la densité g de Y définit elle-même une distribution de Cox. L'une des méthodes de tirage possibles dans ce cas consiste à décomposer la distribution donnée en plusieurs distributions exponentielles (cas notamment de la distribution d'Erlang ou des hyperexponentielles exposées ci-après). Diverses conditions sont nécessaires pour qu'une telle réduction soit possible.

### IV.3.5 - Distribution hyperexponentielle

Une telle distribution a une densité de la forme :

$$ (13) \qquad f(x) = \sum_{i=1}^{p} a_i \lambda_i^x e^{-\lambda_i x} \quad \text{pour } x \ge 0 \text{ avec} $$

$$ (14) \qquad \sum_{i=1}^{p} a_i = 1 $$

Le nombre p de termes dans (13) est l'ordre de l'hyperexponentielle. Les distributions hyperexponentielles sont un cas particulier des distributions de Cox. L'algorithme habituel de tirage indépendant de cette distribution consiste à tirer suivant une distribution exponentielle négative d'espérance $1/\lambda_1$ avec une probabilité $a_1$, ou $1/\lambda_2$ avec une probabilité $a_2$, ou $1/\lambda_i$ avec une probabilité $a_i$, etc...
La génération d'un échantillon autocorrélé peut se faire au moyen
de l'algorithme (2). D'après [2, page 452], la distribution
définie par (13) possède bien la propriété d'infinie divisibilité indiquée
plus haut. Nous allons développer plus longuement ci-après le cas d'une
hyperexponentielle d'ordre 2. Soit :

$$
(15) \quad f(x) = \beta \lambda e^{-\lambda x} + (1-\beta) \mu e^{-\mu x} \quad \text{avec } 0 < \beta < 1
$$

la densité d'une telle distribution. La transformée de Laplace de f est :

$$
f^*(s) = \frac{\beta\lambda}{\lambda+s} + (1-\beta)\frac{\mu}{\mu+s}
$$

Appliquant (12), il vient :

$$
g^*(s) = \frac{\frac{\beta\lambda}{\lambda+s} + (1-\beta) \frac{\mu}{\mu+s}}{\frac{\beta\lambda}{\lambda+\rho s} + (1-\beta) \frac{\mu}{\mu+\rho s}}
$$

ce qui donne après calculs :

$$
g^*(s) = \left( \frac{\lambda + \rho s}{\lambda + s} \right) \left( \frac{\mu + \rho s}{\mu + s} \right) \left( \frac{\lambda \mu + s[\beta \lambda + (1-\beta)\mu]}{\lambda \mu + \rho s[\beta \lambda + (1-\beta)\mu]} \right)
$$

Posons :

$$
(16) \quad \ell = \frac{\lambda \mu}{\rho [\beta \lambda + (1-\beta) \mu]}
$$

$g^*(s)$ peut se décomposer en éléments simples sous la forme :

$$
g^*(s) = A + \frac{B\lambda}{\lambda+s} + \frac{C\mu}{\mu+s} + \frac{D\ell}{\ell+s}
$$

les expressions de A, B, C et D étant, après calculs :

$$
(17) \quad A = \rho
$$

$$
(18) \quad B = \frac{(1-\rho)}{\rho} \left( \frac{\mu-\rho\lambda}{\mu-\lambda} \right) \left( \frac{\rho\ell-\lambda}{\ell-\lambda} \right)
$$
$$ (19) \qquad C = \frac{(1-\rho)}{\rho} \left( \frac{\lambda-\rho\mu}{\lambda-\mu} \right) \left( \frac{\rho\ell-\mu}{\ell-\mu} \right) $$

$$ (20) \qquad D = \frac{(1-\rho)}{\rho} \left( \frac{\lambda-\rho\ell}{\lambda-\ell} \right) \left( \frac{\mu-\rho\ell}{\mu-\ell} \right) $$

La génération d'un échantillon autocorrélé de densité hyperexponentielle d'ordre 2 définie par (15) s'effectue donc au moyen de l'algorithme (2) avec la distribution de Y suivante :

- Y vaut 0 avec la probabilité $A=\rho$

- Avec une probabilité de $(1-\rho)$, Y suit une distribution hyperexponentielle d'ordre 3 de densité :

$$ h(Y) = \frac{B}{1-\rho} \lambda e^{-\lambda Y} + \frac{C}{1-\rho} \mu e^{-\mu Y} + \frac{D}{1-\rho} \ell e^{-\ell Y} $$

où B, C, D et $\ell$ sont donnés par (18), (19), (20), (16).

Nous sommes ramenés au tirage indépendant suivant une distribution hyperexponentielle d'ordre 3 (voir plus haut le rappel de l'algorithme de tirage). Malheureusement, comme on peut le vérifier sur des exemples, les coefficients B, C et D ne sont pas toujours tous trois positifs. Nous sommes amenés à généraliser la forme des distributions hyperexponentielles et à étudier de nouveaux algorithmes de tirage indépendants pour ces cas.

# V - GENERALISATION DE LA FORME DES DISTRIBUTIONS HYPEREXPONENTIELLES

## V.1 - Présentation du problème

Le tirage d'une hyperexponentielle autocorrélée au paragraphe précédent a fait apparaître l'existence de fonctions de densité de la forme :

$$ (21) \qquad f(x) = \sum_{i=1}^{p} a_i \lambda_i^x e^{-\lambda_i x} \qquad \text{avec } x \ge 0 \text{ et} $$

$$ (22) \qquad \sum_{i=1}^{p} a_i = 1 $$
mais où; contrairement à une densité hyperexponentielle classique, certains $a_i$ peuvent être négatifs. Des raisons théoriques (voir précédemment) permettaient cependant d'affirmer dans les exemples précédents que $f(x)$ était bien une densité.

Une fonction $f$ définie sur $[0, \infty[$ est une fonction de densité si et seulement si :

$$ (23) \qquad \int_0^\infty f(x) dx = 1 \quad \text{et} $$

$$ (24) \qquad f(x) \ge 0 \quad \text{pour tout } x \ge 0 $$

Nous appellerons par définition distribution hyperexponentielle généralisée d'ordre $p$ toute distribution dont la densité est donnée par (21), sans autre contrainte sur les coefficients $a_i$ que les conditions (23) et (24) garantissant que $f$ est bien une fonction de densité.

$$ \text{La condition (23) se réduit à (22) : } \sum_{i=1}^{p} a_i = 1 $$

Lorsque tous les coefficients $a_i$ sont positifs, nous sommes ramenés à une hyperexponentielle classique. Dans les autres cas, nous sommes amenés à expliciter les conditions sur les $a_i$ et à étudier de nouveaux algorithmes de tirage. Nous allons étudier successivement les hyperexponentielles généralisées d'ordres 2 et 3.

## V.2 - Hyperexponentielle généralisée d'ordre 2

### V.2.1 - Conditions sur les paramètres de la distribution

La densité définissant une telle distribution est de la forme :
$$
(25) \quad f(x) = a\lambda e^{-\lambda x} + b\mu e^{-\mu x} \quad \text{avec } x \in \mathbb{R}^{+} \text{ et}
$$

$$
(26) \quad a + b = 1
$$

Nous supposons, par un simple choix de notations :

$$
(27) \quad \lambda < \mu
$$

Nous allons rechercher les conditions sur a et b pour que f soit
une fonction de densité.

$f(x)$ positif au voisinage de l'infini impose, en raison de (27).

$$
(28) \quad a > 0
$$

Nous limitons désormais notre étude au cas $b < 0$, sinon nous sommes
dans le cas d'une hyperexponentielle classique.

$$
\text{Comme } f(x) = e^{-\mu x} [a\lambda e^{(\mu-\lambda)x} + b\mu]
$$

$f(x)$ est du signe de la fonction :

$$
\varphi(x) = a\lambda e^{(\mu-\lambda)x} + b\mu
$$

Comme $\varphi$ est une fonction monotone et croissante (car $a$ et $(\mu-\lambda)$ sont positifs d'après (27) et (28)), $\varphi$ donc $f$ sont positives pour tout $x \ge 0$ si et seulement si $\varphi(0) \ge 0$ soit

$$
(29) \quad a\lambda + b\mu \ge 0
$$

condition nécessaire et suffisante (avec $a > 0$ si $\lambda < \mu$, ou $b > 0$ si $\lambda > \mu$
d'après (27) et (28)) pour que f soit une fonction de densité.
### V.2.2 - Etude du cas limite

Ce cas est défini par la condition :

$$ (30) \qquad a\lambda + b\mu = 0 $$

Les équations (26) et (30) permettent alors de calculer $a$ et $b$ soit :

$$ a = \frac{\mu}{\mu-\lambda} \text{ et} $$

$$ b = \frac{\lambda}{\mu-\lambda} $$

Soit $G(x)$ la fonction de densité dans ce cas, c'est-à-dire :

$$ (31) \qquad G(x) = \frac{\lambda\mu}{\mu-\lambda} (e^{-\lambda x} - e^{-\mu x}) $$

dont la transformée de Laplace est (après calculs) :

$$ G^*(s) = \frac{\lambda}{\lambda+s} \cdot \frac{\mu}{\mu+s} $$

La densité $G$ est donc le produit de convolution des densités exponentielles négatives $\lambda e^{-\lambda x}$ et $\mu e^{-\mu x}$, d'où une méthode de tirages aléatoires indépendants suivant $G$: $G$ est la densité de la somme de deux variables aléatoires de densités respectives $\lambda e^{-\lambda x}$ et $\mu e^{-\mu x}$.

### V.2.3 - Construction d'un algorithme de tirages dans le cas général à partir du cas limite

Dans le cas où $a$ et $b$ sont tous deux positifs dans (25), l'algorithme de tirage habituel s'applique : choisir avec des probabilités respectives $a$ et $b$ entre l'une des distributions exponentielles $\lambda e^{-\lambda x}$ ou $\mu e^{-\mu x}$. Nous allons construire un algorithme de tirage pour $b < 0$ (en supposant toujours $\lambda < \mu$).
Pour cela, nous allons chercher une décomposition de $f(x)$ sous la
forme :

$$
(32) \quad f(x) = A \lambda e^{-\lambda x} + B \quad \mathcal{G}(x)
$$

où $\mathcal{G}(x)$ est défini par (31).

Par identification des coefficients de $e^{-\lambda x}$ et $e^{-\mu x}$ entre (25) et (32)
et en utilisant (31), il vient :

$$
(33) \quad A = a + \frac{b\mu}{\lambda}
$$

$$
(34) \quad B = \frac{-b(\mu - \lambda)}{\lambda}
$$

D'après (28) et notre hypothèse $b < 0$, nous avons

$$
A > 0 \text{ et } B > 0.
$$

A et B sont d'autre part liés par la relation :

$$
A + B = a + b = 1
$$

Nous en déduisons un algorithme de tirage suivant la densité f dans
ce cas :

- Tirer avec une probabilité A définie par (33) suivant la distribution exponentielle négative $\lambda e^{-\lambda x}$

- Avec une probabilité de (1-A) tirer suivant la distribution définie par (31) c'est-à-dire faire la somme de deux tirages aléatoires indépendants dont les densités respectives sont $\lambda e^{-\lambda x}$ ou $\mu e^{-\mu x}$.
Algorithme équivalent :

Faire la somme de deux tirages aléatoires indépendants X et Y effectués comme suit :

- X a pour densité $\lambda e^{-\lambda x}$

- Y vaut 0 avec une probabilité de A et est tiré avec une probabilité de (1-A) suivant la densité $\mu e^{-\mu x}$.

L'algorithme précédent, ainsi que l'algorithme classique si a et b sont positifs, recouvrent tous les cas possibles.

### V.3 - Hyperexponentielle généralisée d'ordre 3

La densité d'une telle distribution s'écrit :

$$ (35) \qquad f(x) = a\lambda_1 e^{-\lambda_1 x} + b\lambda_2 e^{-\lambda_2 x} + c\lambda_3 e^{-\lambda_3 x} \qquad x \in \mathbb{R}^+ $$

où a, b et c ont des signes quelconques, avec la relation

$$ (36) \qquad a + b + c = 1 $$

Nous supposons, par un simple choix de notations :

$$ (37) \qquad \lambda_1 < \lambda_2 < \lambda_3 $$

Comme $\int_0^\infty f(x)dx = 1$ d'après (24), f est une fonction de densité si et seulement si $f(x) \ge 0$ pour tout $x \ge 0$. En particulier nous avons les conditions nécessaires : $f(0) \ge 0$ et $f(x) \ge 0$ si x tend vers l'infini. Ces deux conditions s'écrivent respectivement :
$$ (38) \quad a\lambda_1 + b\lambda_2 + c\lambda_3 \geq 0 $$

$$ (39) \quad a > 0 $$

f(x) est du signe de

$$ \varphi(x) = a\lambda_1 e^{\frac{(\lambda_3 - \lambda_1)x}{(\lambda_3 - \lambda_1)(\lambda_3 - \lambda_2)x}} + b\lambda_2 e^{\frac{(\lambda_3 - \lambda_2)x}{(\lambda_3 - \lambda_2)(\lambda_3 - \lambda_1)x}} + c\lambda_3 e^{\frac{-(\lambda_3 - \lambda_1)x}{-(\lambda_3 - \lambda_1)(\lambda_3 - \lambda_2)x}} \\ \text{car } f(x) = \varphi(x)e^{\frac{-(\lambda_3 - \lambda_2)x}{-(\lambda_3 - \lambda_2)(\lambda_3 - \lambda_1)x}}. $$

La dérivée de $\varphi(x)$ vaut :

$$ \varphi'(x) = a\lambda_1(\lambda_3 - \lambda_1)e^{\frac{(\lambda_3 - \lambda_1)x}{(\lambda_3 - \lambda_1)(\lambda_3 - \lambda_2)x}} + b\lambda_2(\lambda_3 - \lambda_2)e^{\frac{(\lambda_3 - \lambda_2)x}{(\lambda_3 - \lambda_2)(\lambda_3 - \lambda_1)x}} $$

La fonction $\varphi$ est donc croissante pour $x$ tel que :

$$ (40) \quad e^{\frac{(\lambda_2 - \lambda_1)x}{(\lambda_2 - \lambda_1)(\lambda_2 - \lambda_3)}} \ge \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} $$

D'autre part, $\varphi(+\infty) = +\infty$.

La condition (40) est vérifiée automatiquement pour tout $x \geq 0$ si $b > 0$ (car le second terme de (40) est alors négatif) ou si :

$$ \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} \le 1 \quad \text{soit} \quad a\lambda_1(\lambda_3 - \lambda_1) + b\lambda_2(\lambda_3 - \lambda_2) \ge 0. $$

Sous ces conditions, $\varphi$ est toujours croissante et la condition (38) est suffisante pour que $\varphi$, donc $f$, soit positif pour tout $x \ge 0$. Dans le cas contraire, il nous faut écrire que le minimum de la fonction $\varphi$ est positif ou nul. Ce minimum $m$ est obtenu pour :

$$ x = \frac{1}{\lambda_2 - \lambda_1} \log \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} \quad \text{et vaut} $$
$$m = a\lambda_1(\lambda_3 - \lambda_1) \left[ \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} \right]^{\frac{\lambda_3 - \lambda_1}{\lambda_2 - \lambda_1}} + b\lambda_2(\lambda_3 - \lambda_2) \left[ \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} \right]^{\frac{\lambda_3 - \lambda_2}{\lambda_2 - \lambda_1}} + c\lambda_3$$

Pour récapituler, la fonction f définie par (35) avec la convention (37) est une fonction de densité si et seulement si :

$$ (39) \quad a > 0 \quad \text{et} $$

$$ (38) \quad a\lambda_1 + b\lambda_2 + c\lambda_3 \ge 0 \quad \text{et} \\ \qquad b > 0 $$

$$ (41) \quad \begin{array}{l} \text{ou } a\lambda_1(\lambda_3 - \lambda_1) + b\lambda_2(\lambda_3 - \lambda_2) \ge 0 \\ \text{ou } a\lambda_1(\lambda_3 - \lambda_1) \left[ \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} \right]^{\frac{\lambda_3 - \lambda_1}{\lambda_2 - \lambda_1}} + b\lambda_2(\lambda_3 - \lambda_2) \left[ \frac{b\lambda_2(\lambda_2 - \lambda_3)}{a\lambda_1(\lambda_3 - \lambda_1)} \right]^{\frac{\lambda_3 - \lambda_2}{\lambda_2 - \lambda_1}} + c\lambda_3 \ge 0 \end{array} $$

Dans le cas limite de (38), c'est-à-dire :

$$ (42) \quad f(0) = a\lambda_1 + b\lambda_2 + c\lambda_3 = 0 $$

ce qui entraîne $\varphi(0) = 0$, la condition (41) se simplifie.

Il suffit de remarquer dans ce cas que $\varphi(x) \ge 0$ pour tout $x \ge 0$ si et seulement si $\varphi$ est une fonction croissante pour tout $x \ge 0$, c'est-à-dire si $\varphi'(0) \ge 0$.

D'où la condition :

$$ (43) \quad a\lambda_1^2 + b\lambda_2^2 + c\lambda_3^2 \le 0 $$

Dans le cas où (42) est vérifiée et où d'autre part :

$$ (44) \quad a\lambda_1^2 + b\lambda_2^2 + c\lambda_3^2 = 0 $$
les équations (36), (42) et (44) permettent le calcul de a, b et c en fonction
de $\lambda_1, \lambda_2, \lambda_3$. La transformée de Laplace de f est, après calculs :

$$
f^*(s) = \frac{1}{\lambda_1 + s} \cdot \frac{2}{\lambda_2 + s} \cdot \frac{3}{\lambda_3 + s}
$$

ce qui montre que f est le produit de convolution des trois densités $\lambda_1 e^{-\lambda_1 x}$,
$\lambda_2 e^{-\lambda_2 x}$, $\lambda_3 e^{-\lambda_3 x}$.

Des algorithmes de tirages indépendants suivant la densité f peuvent
être donnés dans les cas pratiques suivants :

- a, b et c sont tous trois positifs. Nous sommes placés dans le cas d'une hyperexponentielle classique.

- L'une des séries de conditions :

$$
\begin{cases} 0 < a < 1 \\ b\lambda_2 + c\lambda_3 \ge 0 \end{cases} \quad \text{ou} \quad \begin{cases} 0 < b < 1 \\ a\lambda_1 + c\lambda_3 \ge 0 \end{cases} \quad \text{ou} \quad \begin{cases} 0 < c < 1 \\ a\lambda_1 + b\lambda_2 \ge 0 \end{cases}
$$

$$
\begin{cases} 0 < b < 1 \\ a\lambda_1 + c\lambda_3 \ge 0 \end{cases} \quad \text{ou} \quad \begin{cases} 0 < c < 1 \\ a\lambda_1 + b\lambda_2 \ge 0 \end{cases}
$$

$$
\begin{cases} 0 < c < 1 \\ a\lambda_1 + b\lambda_2 \ge 0 \end{cases}
$$

est vérifiée. Etudions par exemple le premier cas. f(x) prend alors la forme :

$$
f(x) = a \lambda_1 e^{-\lambda_1 x} + (1-a) \left[ \frac{b}{1-a} \lambda_2 e^{-\lambda_2 x} + \frac{c}{1-a} \lambda_3 e^{-\lambda_3 x} \right]
$$

avec par hypothèse

$$
\frac{b}{1-a} \lambda_2 + \frac{c}{1-a} \lambda_3 \geq 0
$$

La distribution f est le mélange avec des probabilités respectives
de a et (1-a) d'une distribution exponentielle négative de densité
$\lambda_1 e^{-\lambda_1 x}$ et d'une hyperexponentielle généralisée d'ordre 2 de densité

$$
\frac{b}{1-a} \lambda_2 e^{-\lambda_2 x} + \frac{c}{1-a} \lambda_3 e^{-\lambda_3 x} .
$$
# VI - CONCLUSION

Nous venons d'étudier un algorithme de génération d'échantillons autocorrélés de moyenne et de variance données ou de distribution donnée. Cet algorithme est utilisable pour générer des suites autocorrélées positivement et dont l'autocorrélation est à décroissance géométrique. Le cas de l'hyper-exponentielle autocorrélée nous a amenés à définir une nouvelle technique de tirage des combinaisons de lois exponentielles à coefficients éventuellement négatifs. Ces algorithmes de tirage sont d'un intérêt pratique important car il permettent d'utiliser les différents schémas de génération de variables aléatoires dépendantes proposées par Jacobs et Lewis [3]. Ils justifient la poursuite d'études sur ces schémas de génération dont jusqu'à présent, l'intérêt était limité par les difficultés d'utilisation pratique.
REFERENCES

[1] Cox D.A. - "A use of complex probabilities in the theory of stochastic processes" - Proc. Cambridge - Phil. Soc. 51 - (1955) - pp. 313-319.

[2] Feller W. - "An introduction to probability theory and its applications" - 3rd ed. - Vol. II - Wiley - New York - (1968).

[3] Jacobs P.A., Lewis P.A.W. - "A mixed autoregressive moving average exponential sequence and point process (EARMA 1, 1)" -

[4] Knuth D.E. - "The art of computer programming - Seminumerical algorithms - Random numbers" - Vol. 2 - Chapter 3 - Addison Wesley Publishing Company.

[5] Leroudier J. - "Systèmes adaptatifs à mémoire virtuelle" - Thèse de Doctorat d'Etat - Université de Grenoble - (23 juillet 1977).

[6] Smith C.S. - "Multiplicative pseudo-random generators with prime modulus" - Journal of ACM - Vol. 18 - N° 4 - (October 1977).

[7] Leroudier J., Parent M. - "Discret event simulation modelling of computer systems for performance evaluation" - Rapport IRIA n° 177 - (Juin 1976).