M1 MAPI3

Simulations aléatoires

2023-2024
2
Chapitre I

Simulation d’une variable aléatoire-

Méthode de Monte Carlo (2 semaines)

Pré-requis : variables aléatoires, loi des grands nombres, TCL

Le but de ce chapitre est de présenter des méthodes permettant à un programme
informatique de donner des nombres se comportant comme des variables aléatoires qui
suivent une loi de probabilité donnée.
On utilisera ces nombres aléatoires pour divers types de tâches : principalement,
— se faire une idée qualitative du comportement d’un processus aléatoire (via une
représentation graphique, en général) : observer une seule ou quelques simulations
pour voir l’allure « typique » d’une réalisation (par exemple, pour suggérer une
convergence presque sûre) ;
— estimer une quantité (non-aléatoire) qui peut s’écrire comme une probabilité, une
espérance (voire une loi) et peut donc être approchée par des nombres aléatoires
(méthode de Monte-Carlo) ;
— simuler des données fictives pour ensuite éprouver sur celles-ci l’efficacité de mé-
thodes statistiques, ou pour étalonner des tests.

1 Simulations de variables aléatoires uniformes sur [0, 1]

Une remarque génerale est que les nombres générés par ordinateur sans donnée exté-
rieure ne sont pas aléatoires, puisqu’ils correspondent au résultat d’un programme déter-
ministe. Le but est que ces nombres semblent aléatoires, c’est à dire qu’ils satisfont un
certain nombre de tests statistiques que des nombres vraiment aléatoires doivent satisfaire
avec forte probabilité. On parle de nombres pseudo-aléatoires.
Ces suites sont définies en général par récurrence (ou sont déduites d’une suite définie
par récurrence), tel l’exemple historique (et simpliste) des générateurs par congruence

3
4 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

linéaire comme
Xn+1 = 16807Xn ( mod 231 − 1)

(Ici, 2−31 Xn suit approximativement la loi uniforme sur [0, 1]). Le premier terme X0 de la
suite est appelé la graine (en anglais, seed) de l‘algorithme.
Remarque : Il existe également des moyens d’introduire du “vrai hasard” dans la géné-
ration de nombres pseudo-aléatoires, an utilisant des éléments extérieurs au programme :
par exemple, on peut utiliser l’horloge de l’ordinateur, ou des phénomènes physiques (ra-
dioactivité, phénomènes quantiques, etc.)
En général, la graine est fixée en fonction du nombre de millisecondes de l’horloge de
l’ordinateur au moment où l’on commence la simulation (seed(a = None)). Attention, on
n’utilise cette méthode que pour la graine, tandis que les tirages suivants s’en déduiront :
c’est une forme de compromis entre l’aléas pur et le déterminisme.
A partir de maintenant, nous allons supposer que le problème de la génération de
nombres aléatoires est résolu et donc que nous avons le moyen de simuler une suite arbi-
trairement longue de variables aléatoires indépendantes, uniformes sur [0, 1] en utilisant
successivement la même fonction uniform(). On va maintenant décrire le moyen de simuler
des variables suivant d’autres lois.

Exemple : X = (uniform() < p) simule une Bernoulli de paramètre p car X vaut 0 ou

1 avec
P[X = 1] = P[uniform < p] = p.

Simulation de variables binomiales : On utilise le fait qu’une loi binomiale de pa-

ramètre p ∈ [0, 1] et n ≥ 1 est une somme de variables de Bernoulli indépendantes de
paramètre p.

Ainsi, pour simuler une variable aléatoire X ∼ B(n, p), on commence par tirer n va-
riables aléatoires U1 , . . . , Un indépendantes uniformes sur [0, 1], puis on construit les n
variables de Bernoulli Yk = 1[0,p] (Uk ) , pour 1 ≤ k ≤ n et on pose enfin X = Y1 + · · · + Yn .

2 Simulations de variables aléatoires sur R

2.1) Variables aléatoires discrètes

Proposition 2.1. Soit x1 , x2 , · · · , xn des nombres réels tous différents et soit p1 , p2 , · · · , pn
des nombres réels positifs tels que ni=1 pi = 1. On pose s0 = 0 et pour tout 1 ≤ k ≤ n,
P
2. SIMULATIONS DE VARIABLES ALÉATOIRES SUR R 5
Pk
sk = i=1 pi . Soit U est une variable aléatoire de loi uniforme U([0, 1]) et
n
X
X= xk 1(sk−1 ≤U ≤sk ) .
k=1

Alors, X est une variable aléatoire de loi discrète P = p1 δx1 + p2 δx2 + · · · + pn δxn .

Démonstration. On vérifie que pour tout i ∈ {1, . . . , n},

P(X = xi ) = P(si−1 ≤ U ≤ si ) = si − si−1 = pi .

Variables aléatoires discrètes

Si on veut simuler la loi p1 δx1 + p2 δx2 + · · · + pn δxn en utilisant le vecteur des pi et
celui des xi
(si )i = sommes cumulees des (pi )i ;
U = uniform();
i = 1;
tant que U > si , i = i + 1;
sortir xi ;

2.2) Calcul de lois

Si on X a la même loi que f (U1 , . . . , Un ) avec U1 , . . . , Un indépendantes et uniformes,
alors on peut simuler X :
Simuler X de loi f (U1 , . . . , Un )
On applique f à des variables aléatoires uniformes U1 , . . . , Un

Exemple 1 : v.a. uniformes sur un intervalle quelconque

Proposition 2.2. Si U est une variable uniforme sur [0, 1], alors a+(b−a)U est uniforme
sur [a, b].

Démonstration. Il suffit de remarquer que

P(a + (b − a)U ≤ t) = P(U ≤ (t − a)/(b − a)) = (t − a)/(b − a).

Exemple 2 : On veut simuler une v.a. à valeurs dans [0, 1] et de densité 2x. On se rend
compte que si U1 et U2 sont deux variables indépendantes et uniformes sur [0, 1], alors
max(U1 , U2 ) a la bonne densité (calculer la fonction de répartition). On en déduit qu’il
suffit de prendre le maximum entre deux uniformes pour simuler notre v.a.
6 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

Changement de variable : De manière générale, pour calculer la loi de f (U1 , . . . , Un ),

qui est une mesure image, il est utile de rappeler la formule du changement de variable
Z b Z ϕ(b)
0
f (ϕ(t))ϕ (t)dt = f (x)dx,
a ϕ(a)

qu’on utilise formellement en disant qu’on effectue le changement de variable x = ϕ(t) et

que par conséquent dx = d(ϕ(t)) = ϕ0 (t)dt.

Exemple 3 : Si V est uniforme sur [0, π/2], calculons la loi de sin2 (V ). Pour toute
fonction mesurable bornée g, on a

2 π/2
Z
2
E[sin (V )] = g(sin2 (θ))dθ
π 0
2 1 √

Z
= g(x)d arcsin( x)
π 0
2 1
Z
1 1
= g(x) √ q dx
π 0 2 x √ 2
1− x
Z 1
1
= g(x) p dx.
0 π x(1 − x)

On en déduit que si U est uniforme sur [0, 1], sin2 (πU/2) suit la loi sur [0, 1], de densité
p
x 7→ 1/(π x(1 − x)), c’est à dire la loi beta de paramètres 1/2 et 1/2. Cette loi est
également appelée loi de l’arc sinus, car sa fonction de répartition fait intervenir la fonction
x 7→ arcsin x.

Exemple 4 : Simulation de variables aléatoires uniformes sur le disque unité. En coor-

données cartésienne, la mesure uniforme sur le disque unité peut être donnée par l’expres-
sion (1/π)1x2 +y2 ≤1 dxdy, soit en coordonnées polaires (1/π)1r≤1 rdrdθ. En posant s = r2 ,
on obtient (1/2π)1s≤1 dsdθ puisque ds = 2rdr.
Ainsi, pour (x, y) uniforme sur le disque unité, s = r2 et θ sont indépendant, uniformes
respectivement sur [0, 1] et [0, 2π).

Proposition 2.3. Si U et V sont deux variables indépendantes uniformes sur [0, 1], alors
√ √
( U cos(2πV ), U sin(2πV )) est uniforme sur le disque unité.

2.3) Méthode d’inversion de la fonction de répartition

Soit µ une loi de probablité sur R (muni de la tribu borélienne), et F sa fonction de
répartition :
F (x) := µ((−∞, x]).
2. SIMULATIONS DE VARIABLES ALÉATOIRES SUR R 7

On rappelle que F est croissante, continue à droite, admet une limite à gauche en tout
point, tend vers 0 en −∞ et vers 1 en +∞.

Définition 2.4. Pour t ∈ (0, 1), on définit la pseudo-inverse G de F (appelée aussi

fonction de quantile de µ) par

G(t) := inf{x ∈ R, F (x) ≥ t} < ∞.

On remarque en particulier que G(t) est finie à cause des limites de F en ±∞.
Quelques propriétés

Proposition 2.5. 1. G(t) ≤ x ⇔ F (x) ≥ t

2. Si F restreint à ]a, b[ établit une bijection de ]a, b[ vers ]0, 1[, (d’inverse F −1 ), alors
G = F −1 sur ]0, 1[.

Démonstration. 1. Par définition de G, si on a F (x) ≥ t alors G(t) ≤ x. Réciproque-

ment, si G(t) ≤ x, alors pour tout ε > 0, on a G(t) < x + ε, ce qui implique que
t ≤ F (x + ε). Comme F est continue à droite, on peut faire tendre ε vers 0 et
obtenir t ≤ F (x).
2. Pour t ∈]0, 1[, F −1 (t) est l’unique réel x tel que F (x) = t, et par croissance, c’est
donc le plus petit réel x tel que F (x) ≥ t. On a donc F −1 (t) = G(t) par définition
de G.

Proposition 2.6. Soit U une variable aléatoire uniforme sur [0, 1], alors, la variable
aléatoire G(U ) suit la loi µ.

Démonstration. Si U est uniforme sur [0, 1], U est presque sûrement dans (0, 1), donc la
variable aléatoire G(U ) est bien définie.
Calculons maintenant la fonction de répartition de la variable G(U ).

P[G(U ) ≤ x] = P[U ≤ F (x)] = F (x).

Donc G(X) suit la loi µ.

Inversion de la fonction de répartition

Soit maintenant une distribution sur R
— On calcule sa fonction de répartition F ,
— On calcule la fonction G
— On applique G à des variables aléatoires uniformes G(Uk )
8 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

Corollaire 2.7. Si (Xn )n≥1 sont indépendantes, uniformes sur [0, 1], alors (G(Xn ))n≥1
sont indépendantes uniformes de loi µ.

Démonstration. En effet, si l’on applique la même fonction mesurable à une suite de

variables aléatoires i.i.d. on obtient une suite de variables aléatoires i.i.d.
Comme la fonction G est mesurable, le résultat est immédiat.

La méthode présentée ici s’applique à n’importe quelle distribution à valeurs dans R,

cependant elle demande d’inverser la fonction de répartition, ce qui peut être malaisé (par
exemple, pour une variable gaussienne). Nous allons voir certains cas particuliers où cette
méthode peut être utilisée.

Exemple : Simulation de variables exponentielles

Proposition 2.8. Si Un est une variable uniforme sur [0, 1], et λ > 0, alors on a
−(1/λ) log(Un ) est une suite de v.a. indépendantes de loi exponentielle de paramètre λ.

Démonstration. On rappelle que la fonction de répartition d’une variable exponentielle

est
F (t) = 1 − e−λt , t ≥ 0.

Ainsi, si on inverse la relation on obtient

1
u = 1 − e−λt ⇔ t = − log(1 − u).
λ
On en déduit donc que si U suit une loi uniforme alors −1/λ log(1 − U ) suit une loi
exponentielle.
Il reste à remarquer que si U suit une loi uniforme sur [0, 1] alors 1 − U également, on
en déduit donc la proposition.

Alternativement on vérifie que pour tout x ≥ 0,

P[−(1/λ) log(Xn ) ≤ x] = P[Xn ≥ e−λx ] = 1 − e−λx ,

ce qui prouve que −(1/λ) log(Xn ) est une variable exponentielle de paramètre λ.

2.4) Mélange de lois

Soit N, X1 , X2 trois v.a. réelles indépendantes. On suppose que N est à valeurs dans
{1, 2}, et que X1 et X2 sont à densité. Calculons la loi de XN .
2. SIMULATIONS DE VARIABLES ALÉATOIRES SUR R 9

On note f1 et f2 les densités de X1 et X2 . Pour toute fonction mesurable bornée g,

comme g(XN ) = 11N =1 g(X1 ) + 11N =2 g(X2 ), on calcule

E[g(XN )] = E[11N =1 g(X1 ) + 11N =2 g(X2 )]

= E[11N =1 ]E[g(X1 )] + E[11N =2 ]E[g(X2 )]
Z Z
= P[N = 1] g(x)f1 (x)dx + P[N = 2] g(x)f2 (x)dx
Z R R

= g(x) (α1 f1 (x) + α2 f2 (x)) dx

R

ce qui montre que XN a pour densité x 7→ α1 f1 (x) + α2 f2 (x). Plus généralement, on a

Proposition 2.9. Si N, X1 , · · · , Xn sont des v.a. réelles avec N à valeurs dans {1, · · · , n}
de masses α1 , . . . , αn ≥ 0 et indépendante des autres, alors la v.a. XN a pour loi

α1 PX1 + . . . + αn PXn .

Pour simuler une variable aléatoire de loi µ, où µ = α1 µ1 +· · ·+αn µn , avec α1 , . . . , αn ≥

0, α1 + · · · + αn = 1, on peut donc faire comme ceci :

Mélange de loi
On simule une réalisation k de N à valeurs dans {1, · · · , n} de probas α1 , . . . , αn ,
puis on simule X de loi µk .

Par exemple, la loi µ sur R de densité ρ par rapport à la mesure de Lebesgue, avec
2 √ 2 √
ρ(x) = e−x /2 / 8π si x ∈/ [0, 1] et ρ(x) = (1/2) + e−x /2 / 8π si x ∈ [0, 1], est la demi-
somme de la loi uniforme sur [0, 1] et de la loi gaussienne centrée réduite. On en déduit
que si Z1 , Z2 , J sont des variables indépendantes, respectivement uniforme sur [0, 1],
gaussienne centrée réduite et uniforme sur {1, 2}, alors ZJ suit la loi µ.
On peut généraliser la proposition précédente pour un mélange infini de lois :

Proposition 2.10. Si (Xt )t∈R sont des v.a. réelles de densités ft , et si T est une variable
v.a. réelle de densité g, et indépendante des (Xt )t∈R , alors la v.a. XT a pour densité
Z
x 7→ ft (x)g(t)dt.
R

Pour simuler une variable aléatoire de densité

Z
x 7→ ft (x)g(t)dt,
R

on peut donc faire comme ceci :

Mélange de loi
On simule une réalisation t de T , puis on simule X de densité ft .
10 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

2.5) Méthode de rejet

La méthode du rejet permet de simuler une variable aléatoire de loi η en conditionnant
une loi µ que l’on sait simuler facilement.

a) Conditionnement par un évènement

Pour obtenir une simulation de Z sachant un événement A qui dépend de Z (on

dit que c’est un évènement σ(Z)-mesurable), il suffit de simuler de manière répétée et
indépendante Z et de rejeter les résultats tant que A n’est pas réalisé (Z peut être une
variable aléatoire multi-dimensionnelle).
Proposition 2.11. Soit Z une v.a. et A un évènement de probabilité non nul. Considérons
(Zn , An )n≥1 une suite d’éléments aléatoires indépendants de même loi que (Z, A). Notons
K = inf{n ≥ 1 : An est réalisé} : alors la v.a. ZK a pour loi la loi conditionnelle de Z
sachant A.
Démonstration. Pour tout borélien B, on a
X
P[ZK ∈ B] = P[Zn ∈ B; Ac1 ; . . . ; Acn−1 ; An ]
n≥1
X P[Z ∈ B; A]
= (1 − P[A])n−1 · P[Zn ∈ B; An ] = = P[Z ∈ B|A]
n≥1
P[A]

Conditionnement par A
On répète la même simulation et on garde uniquement les valeurs qui vérifient A
On a
P[K = n] = (1 − P[A])n−1 · P[A],
donc K est un temps aléatoire qui suit une loi géométrique de paramètre P[A] (espérance
1/P[A]). Donc plus A est probable, et plus la simulation est rapide.

Exemple : Par définition, la loi uniforme sur un ensemble de mesurable B de Rd de

1
mesure de Lebesgue 0 < |B| < +∞ est la mesure de densité |B| 1B . Pour tout A ∈ Rd , la
mesure de A sous la loi uniforme sur B est donc donnée par
|A ∩ B|
Z
1
1B (x)dx = .
A |B| |B|
Maintenant, prenons une variable aléatoire X uniforme dans B, et calculons la loi condi-
tionnelle de X sachant que X ∈ C ⊂ B. Pour tout borélien A,
 −1
P[X ∈ A ∩ C] |A ∩ C ∩ B| |C ∩ B| |A ∩ C|
P[X ∈ A|X ∈ C] = = = ,
P[X ∈ C] |B| |B| |C|
2. SIMULATIONS DE VARIABLES ALÉATOIRES SUR R 11

donc la loi de X uniforme sur B sachant que X ∈ C ⊂ B est en fait la loi uniforme sur
C.

Simulation de variables uniformes sur le disque unité : Nous allons utiliser la

méthode de rejet, en observant que la loi uniforme sur le disque unité est la loi uniforme
sur le carré [−1, 1]2 , conditionnée par le fait que la variable prend ses valeurs sur le disque.
La loi uniforme sur [−1, 1]2 est aisément simulée en prenant deux variables indépendantes
(X, Y ), uniformes sur [−1, 1]. Dans cet exemple, la loi auxilliaire µ s’écrit

µ(dx, dy) = 1[0,1] (x)1[0,1] (y)dxdy

et η = µ(·|A) où
A = (x2 + y 2 ≤ 1 .


On rejette donc avec probabilité 1 − π/4 ' 0.21 soit un peu plus de 20 pourcent des
indices.

Simulation de variables uniformes sur la boule unité et sur la sphère unité :

La méthode de rejet s’applique pour la mesure uniforme sur la boule unité en dimension
N pour tout N ≥ 1, à partir de la mesure uniforme sur l’hypercube [−1, 1]N , obtenue en
prenant N variables indépendantes uniformes sur [−1, 1].
Il y a cependant un problème pour N grand car la probabilité de garder un indice
k tend rapidement vers 0 quand N tend vers l’infini. Cette probabilité est le quotient
du volume de la boule par celui de l’hypercube, et on peut montrer que cela donne
π N/2 /2N Γ(N/2 − 1), qui tend vers zéro plus qu’exponentiellement vite.
Enfin, si on a des variables uniformes sur la boule unité de RN , on obtient des uniformes
sur la sphère unité en divisant les uniformes sur la boule par leur norme (admis ici).

b) Méthode du rejet

Ici, nous supposons que la loi de la variable aléatoire X d’intérêt (éventuellement

multidimensionnelle) possède une densité f connue, mais dont la simulation directe n’est
pas facile. Le principe de la méthode consiste à simuler une autre variable aléatoire Y de
densité g et de conditionner par un événement A tel que la loi de Y conditionnelle est celle
de X. Cette idée remonte à Von Neumann en 1947. La propriété s’énonce précisément
ainsi.

Proposition 2.12. Soit X et Y deux v.a. de densités respectives f et g telles que

f (x) ≤ c · g(x)
12 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

pour une certaine constante c.

Soit U une variable aléatoire de loi uniforme sur [0, 1]. Alors la loi de Y sachant
{cU g(Y ) < f (Y )} est la loi de X. De plus,

P[cU g(Y ) < f (Y )] = 1/c.

Démonstration. Posons A = {cU g(Y ) < f (Y )}. Pour tout borélien B, on a

P[Y ∈ B|A] = P[Y ∈ B; A]/P[A]

Z
1
= g(y)dydu
P[A] (y,u)∈Rd ×[0,1]:y∈B,c u g(y)<f (y),g(y)6=0
Z
1 f (y)
= g(y) dy
P[A] y∈Rd :y∈B,g(y)6=0 c g(y)
Z
1
= f (y)dy.
c P[A] y∈B

Comme P[A] = 1/c (en prenant B = Rd dans l’égalité), on trouve bien

Z
P[Y ∈ B|A] = f (y)dy.
y∈B

Méthode du rejet
On simule plusieurs fois Y de densité g et U uniforme jusqu’à ce que cU g(Y ) ≤ f (Y ).
La dernière variable Y suit la loi de X de densité f .

Pour avoir un algorithme aussi rapide que possible, il est souhaitable de ne pas trop rejeter
d’indices, et donc que c soit aussi petit que possible. L’optimum est en prenant pour c le
maximum de ρ mais dans certains cas il n’est pas facile à calculer, on se contente alors
d’un majorant “raisonnable”.

Simulation de variables aléatoires gaussiennes : Pour simuler une variable gaus-

sienne centrée réduite, il peut être malaisé d’utiliser la fonction de répartition car elle
n’a pas d’expression en termes de fonctions “usuelles”. On va donner deux méthodes de
simulations, d’abord par la méthode du rejet, puis en utilisant un changement de variables.

Pour la méthode de rejet on peut choisir comme loi auxiliaire, la loi d’un produit
d’une variable exponentielle de paramètre λ > 0 (simulée précédemment) par une variable
uniforme sur {−1, 1}. Cette loi a pour densité par rapport à la mesure de Lebesgue
λ −λ|x|
f (x) = e .
2
2. SIMULATIONS DE VARIABLES ALÉATOIRES SUR R 13

En effet,
1 1
E(φ(sE)) = E(φ(−E)) + E(φ(E))
2Z 2
1 ∞ 1 ∞
Z
= φ(−x)λ exp(−λx)dx + φ(x)λ exp(−λx)dx
2 0 2 0
Ainsi, la densité r
f (x) 2 1 λ|x|−x2 /2
ρ(x) = = e ,
g(x) πλ
dont le maximum est obtenu pour x = q±λ 2(étudier le logarithme de ρ). Ce maximum
λ /2
peut être pris pour c et on a alors c = π2 e λ .
Afin d’optimiser la méthode, il faut choisir c pour avoir une probabilité de rejet la
plus faible possible. On rappelle que la probabilité de rejet vaut 1 − 1/c, c’est à dire qu’on
souhaite minimiser c (ou de manière équivalente, on minimise son logarithme). Le choix
de λ optimisant c est λ = 1, qui donne c ' 1.31 et une probabilité de rejeter de l’ordre
de 24%.

On peut également utiliser la méthode suivante, due à Box and Müller.

Proposition 2.13 (Méthode de Box and Müller). Si U, V sont des indépendantes uni-
formes sur [0, 1], alors les variables aléatoires définies par
p p
−2 log(U ) cos(2πV ) et −2 log(U ) sin(2πV )

sont deux variables indépendantes gaussiennes centrées réduites.

Démonstration. Montrons que si (X, Y ) est un couple de variables gaussiennes centrées

réduites indépendantes alors
(loi) p p
(X, Y ) = ( −2 log(U ) cos(2πV ), −2 log(U ) sin(2πV ),

où U et V sont des variables aléatoires indépendantes uniformes sur [0, 1].

Pour cela, nous allons faire des changements de variables successifs dans la loi de (X, Y ) :
2 +y 2 )/2
(1/2π)e−(x dxdy.

En coordonnées polaires, cela donne

2 /2
(1/2π)e−r rdrdθ,

soit pour s = r2 /2,

(1/2π)e−s dsdθ.
14 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

On voit ainsi que l’on obtient ainsi la loi de deux variables indépendantes, l’une S ex-
ponentielle de paramètre 1 et l’autre uniforme θ sur [0, π). Pour conclure il suffit de
rappeller qu’on peut obtenir une variable exponentielle à partir d’une variable uniforme
√
U par − log(U ). En utilisant le changement de variable x = r cos(θ) = 2s cos(θ) et
√
y = 2s sin(θ), on en déduit que si U et V sont deux variables indépendantes, uniforme
p p
sur [0, 1], alors −2 log(U ) cos(2πV ) et −2 log(U ) sin(2πV ) sont deux variables gaus-
siennes centrées réduites indépendantes.

Toutes les variables gaussiennes à valeurs réelles peuvent être immédiatement obtenues
à partir de variables gaussiennes centrées réduites.

Proposition 2.14. Soit X est une variable gaussienne centrée réduite, alors σX + µ est
une variable gaussienne de moyenne µ et de variance σ 2 .

Démonstration. Soit f une fonction continue bornée, alors

Z ∞  2
1 x
E(f (m + σX)) = f (m + σx) √ exp − dx
−∞ 2π 2
Z ∞
(y − m)2
 
1
= f (y) √ exp − dy
−∞ 2πσ 2 2σ 2

3 Méthode de Monte Carlo

3.1) Principe de la méthode
La méthode de Monte-Carlo consiste à obtenir une valeur approchée de l’espérance
E[X] d’une variable aléatoire (intégrable) X en simulant un grand nombre N de variables
aléatoires X1 , . . . , XN indépendantes et de même loi que X, afin d’appliquer la loi forte
des grands nombres :
X1 + . . . + XN
E[X] ' .
N
Ne pas oublier que pour une fonction mesurable ϕ telle que ϕ(X) est intégrable, les
variables aléatoires ϕ(X1 ), . . . , ϕ(XN ) sont indépendantes et de même loi que ϕ(X), donc
la loi forte des grands nombres s’écrit alors
ϕ(X1 ) + . . . + ϕ(XN )
E[ϕ(X)] ' .
N
En particulier, pour un mesurable A, l’égalité précédente appliquée à l’indicatrice ϕ = 1A
nous donne que
1  
E[ϕ(X)] = P[X ∈ A] ' Card {n entre 1 et N tels que Xn ∈ A} .
N
3. MÉTHODE DE MONTE CARLO 15

Ceci rejoint parfaitement l’intuition qu’on a d’une probabilité P[X ∈ A] : c’est la propor-
tion de X qui tombe dans A lorsqu’on tire un grand nombre de fois X.

Méthode de Monte-Carlo
On simule X1 , . . . , XN indépendantes de même loi que X. La proportion
1  
Card {n entre 1 et N tels que Xn ∈ A}
N
d’incides tels que Xn ∈ A est une approximation de P[X ∈ A].

Exemple : Si on connaît la mesure |B| d’un domaine B ⊂ Rn , pour estimer la mesure

|A| d’un domaine A ⊂ B, on peut considérer un point uniforme X dans B et observer
que
P[X ∈ A] = |A|/|B|.

Pour estimer
 |A|, il suffit donc d’estimer P[X
 ∈ A] ce qui revient à calculer la proportion
1
N
Card {n entre 1 et N tels que Xn ∈ A} pour Xn points tirés uniformément dans |B|.

3.2) Intervalle de confiance

Si ϕ(X)2 est intégrable, l’erreur d’estimation s’obtient par le théorème central limite :
√
 
ϕ(X1 ) + . . . + ϕ(XN )
N − E[ϕ(X)]
N

converge en loi vers une gaussienne centrée de variance σ 2 = Var(ϕ(X)). En particulier,

 
ϕ(X1 ) + . . . + ϕ(XN ) aσ
P − E[ϕ(X)] ≤ √ ' P [|Z| ≤ a] ,
N N
où Z est une gaussienne centrée réduite. On en déduit alors l’intervalle de confiance
suivant : avec probabilité P [|Z| ≤ a], la vraie valeur E[ϕ(X)] appartient à l’intervalle
aléatoire de centre ϕ(X1 )+...+ϕ(X
N
N)
et de demi-largeur √aσN .
En pratique, on fixe un niveau de confiance P [|Z| ≤ a] proche de 1 (par exemple 0, 95),
et on cherche le a qui nous donne ce niveau à l’aide de tables de valeurs (par exemple, pour
a = 1, 96, P [|Z| ≤ a] ' 0.95). Ceci nous donne : avec probabilité 0, 95, la vraie valeur
E[ϕ(X)] appartient à l’intervalle aléatoire de centre ϕ(X1 )+...+ϕ(X
N
N)
et de demi-largeur
1,96σ
√ .
N
L’écart-type σ n’est pas forcément connu, et il faut parfois l’estimer, ou au moins le
majorer, pour assurer le niveau de confiance. Dans le cas où on estime une probabilité
P[X ∈ A] (c’est le cas particulier où ϕ = 1A ,), la variance est celle d’une Bernoulli, que
l’on peut toujours majorer par 1/4. On en déduit alors l’intervalle de confiance suivant :
16 CHAPITRE I. SIMULATION DE V.A. ET MÉTHODE DE MC (2 SEMAINES)

avec probabilité supérieureà P [|Z| ≤ a], la vraie valeur P[X ∈A] appartient à l’intervalle
aléatoire de centre N1 Card {n entre 1 et N tels que Xn ∈ A} et de demi-largeur 2√aN .
Par exemple, pour a = 1, 96 de telle sorte que P [|Z| ≤ a] ' 1, 96, ceci donne :

Intervalle de confiance pour Monte-Carlo

On simule X1 , . . . , XN indépendantes de même loi que X. Avec probabilité supé-
rieure à 0, 95, la quantité P[X ∈ A] appartient
 à l’intervalle aléatoire de centre
N
1
Card {n entre 1 et N tels que Xn ∈ A} et de demi-largeur 0,98
√ .
N

Exemple : Pour revenir à l’exemple de la section précédente, où on estime l’aire |A|

d’un domaine A ⊂ B ⊂ Rn , on obtient donc un intervalle de confiance pour

P[X ∈ A] = |A|/|B|
 
0,98
centré en N1 Card {n entre 1 et N tels que Xn ∈ A} et de demi-largeur √ .
N
Pour en
faire un intervalle de confiance pour |A|, il suffit de le multiplier par |B|.
Chapitre II

Modèle gaussien (2 semaine)

1 Définitions et simulations
Définition 1.1. Une variable aléatoire gaussienne centrée réduite est une variable
aléatoire réelle de densité  2
1 x
f (x) = √ exp − .
2π 2
Une variable aléatoire gaussienne est une variable aléatoire X = µ + σZ avec µ, σ ∈ R
et Z une variable aléatoire gaussienne centrée réduite. On note X ∼ N (µ, σ 2 ). Si σ 6= 0,
la densité de X est alors

(x − µ)2
 
1
f (x) = √ exp − .
2πσ 2 2σ 2

Rappelons que pour simuler une variable aléatoire centrée réduite, on peut utiliser
la méthode de Box-Müller : si U, V sont des indépendantes uniformes sur [0, 1], alors les
variables aléatoires définies par
p p
−2 log(U ) cos(2πV ) et −2 log(U ) sin(2πV )

sont deux variables indépendantes gaussiennes centrées réduites. Ensuite, en multipliant

par σ et en ajoutant µ, on obtient une une variable aléatoire gaussienne de moyenne µ et
de variance σ 2 .

Définition 1.2. Un vecteur aléatoire X = (X1 , . . . , Xn ) est un vecteur gaussien si toute

combinaison linéaire de ses composantes est une variable aléatoire gaussienne. La moyenne
µ de X, et la matrice de covariance Σ de X caractérisent la loi de X, et on note X ∼
N (µ, Σ).

17
18 CHAPITRE II. MODÈLE GAUSSIEN (2 SEMAINE)

Pour rappel, le vecteur moyenne µ est

 
E[X1 ]
..
E[X] = 
 
. 
E[Xn ]
et la matrice de covariance Σ est la matrice n × n

Var(X) = E (X − E[X]) · (X − E[X])T

 
  
X1 − E[X1 ]
..
 
= E   · X1 − E[X1 ], · · · , Xn − E[Xn ]  ,
  
.
Xn − E[Xn ]

(autrement dit, Σi,j = cov(Xi , Xj )).

Attention ! Il ne suffit pas que toutes les composantes d’un vecteur soit des gaussiennes
pour que le vecteur soit gaussien. Pour montrer qu’un vecteur est gaussien, on peut mettre
en place plusieurs stratégies :
— on l’exprime comme combinaisons linéaires d’un autre vecteur gaussien,
— ou bien on calcule la densité. Si det(Σ) 6= 0, la densité de X par rapport à la
mesure de Lebesgue sur Rn est alors
(x − µ)T Σ−1 (x − µ)
 
1
f (x) = p exp − ,
(2π)d det(Σ) 2

où (x − µ) et (x − µ)T désignent respectivement un vecteur colonne et un vecteur

ligne.
— Une autre possibilité est de calculer la fonction caractéristique. La fonction carac-
téristique d’un vecteur gaussien X de moyenne µ et de matrice de covariance Σ est
donnée par
ϕX (z) (= E[eihz,Xi ]) = exp(ihz, µi − hz, Σzi),
ce qui est l’outil principal pour montrer des égalités en loi, notamment celles-ci.

Proposition 1.3. 1. Si X ∼ N (µ, Σ) (en colonne), alors AX+b ∼ N (Aµ+b, AΣAT )

2. Si X ∼ N (µ, Σ), alors pour tout 1 ≤ i, j ≤ n, on a Xi et Xj indépendants si et
seulement si Σi,j = 0.

Démonstration. Il suffit à chaque fois de calculer la fonction caractéristique.

On peut donc utiliser la simulation de variables gaussiennes indépendantes pour obtenir

un vecteur X ∼ N (µ, In ). Ensuite, on peut le multiplier par A tel que AAT = Σ pour
obtenir AX ∼ N (µ, AAT ). La décomposition Σ = AAT peut se faire par exemple à
l’ordinateur via la décomposition de Cholesky.
2. ESTIMATEUR ET CONDITIONNEMENT 19

2 Estimateur et conditionnement
Etant donné un couple de v.a. (X, Y ), on souhaite estimer X en observant Y . Une
stratégie naturelle est de prendre, parmi les fonctions f (Y ) de Y , la variable aléatoire
f (Y ) qui minimise l’erreur quadratique

kX − f (Y )kL2 := E[(X − f (Y ))2 ]1/2 .

On note L2 (Y ) l’espace vectoriel {f (Y ) : fonction f telle que f (Y ) soit de carré intégrable}.

On cherche donc Z ∈ L2 (Y ) tel que kX − ZkL2 soit minimal.
La réponse est donnée par la projection orthogonale de X sur L2 (Y ) dans l’espace de
Hilbert des variables de carré intégrable. En effet, en utilisant le produit scalaire

hX, ZiL2 := E[XZ],

on est dans le cadre suivant.

Proposition 2.1 (Projection dans un Hilbert). Soit H, h·, ·i, k · k un espace de Hilbert
et S un sous-espace vectoriel fermé de H. Alors pour tout x ∈ H, il existe un unique
p(x) ∈ S tel que
kx − p(x)k = inf kx − zk.
z∈S

De plus il vérifie ∀z ∈ S, hx − p(x), zi = 0 et cette propriété le caractérise. On dit que

p(x) est la projection orthogonale de x sur S.

Définition 2.2. Etant donné une variable aléatoire X et un vecteur Y de variables aléa-
toires, on appelle espérance conditionnelle de X sachant Y la projection orthogonale de
X sur L2 (Y ), c’est-à-dire la seule variable aléatoire E(X|Y ) ∈ L2 (Y ) qui vérifie


E (X − E(X|Y ))Z = 0

pour tout Z ∈ L2 (Y ).

L’espérance conditionnelle minimise l’erreur quadratique kX − E[X|Y ]kL2 . Bien en-

tendu, si X est un vecteur colonne, on définit
 
E[X1 |Y ]
..
E[X|Y ] =  .
 
.
E[Xn |Y ]

On est souvent amener à considérer des vecteurs gaussiens (X1 , . . . , Xn , Y1 , . . . Ym ) ∈ Rn+m

et à n’en observer qu’une partie Y = (Y1 , . . . , Ym ). Un estimateur naturel de Xi est
alors donnée par l’espérance conditionnelle E[Xi |Y ] de Xi sachant Y . Les calculs sont
particulièrement agréables comme le montre les propriétés suivantes :
20 CHAPITRE II. MODÈLE GAUSSIEN (2 SEMAINE)

Proposition 2.3. Soit (X, Y ) = (X1 , . . . , Xn , Y1 , . . . Ym ) ∈ Rn+m un vecteur gaussien.

Alors la variable aléatoire
! E[Xi |Y ] est une fonction affine de Y = (Y1 , . . . , Ym ).
X
De plus, si ∼ N (µ, Σ) (en vecteurs colonnes) avec
Y
! !
µX ΣX ΣXY
µ= , Σ= ,
µY ΣY X ΣY
on a
E[X|Y ] = µX + ΣXY Σ−1
Y (Y − µY ).

La première partie nous dit que la projection est en fait une projection sur un espace
vectoriel de dimension fini, ce qui rend les calculs plus faciles.

Démonstration. Commençons par vérifier que le vecteur X − µX + ΣXY Σ−1

Y (Y − µY )
est centré. Ainsi sa covariance avec une variable aléatoire Z s’écrit

Cov(X − µX + ΣXY Σ−1 −1



Y (Y − µY ) , Z) = E (X − µX + ΣXY ΣY (Y − µY ) )Z .

On vérifie que les coordonnées de X − µX + ΣXY Σ−1

Y (Y − µY ) sont orthogonales à
Y1 , . . . , Yn en calculant

X − µX + ΣXY Σ−1

T
E Y (Y − µY ) Y

= E(XY T ) − µX E(Y T ) + ΣXY Σ−1 T

Y E((Y − µY )Y )

= ΣXY + µx µTY − µX µTY + ΣXY Σ−1

Y (ΣY )

=0

Grace à la remarque précédente, on en déduit que les deux vecteurs ont une coviariance
nulle, et comme le vecteur est gaussien, cela signifie que les coordonnées sont indépen-
dantes des Yi .
Chaque coordonnée est alors indépendantes de toute variable aléatoire Z ∈ L2 (Y ). On
a donc


E (Xi − E(Xi |Y ))Z = 0
pour tout Z ∈ L2 (Y ), comme voulu.

Plus généralement, on peut montrer que

!
X
Proposition 2.4. Soit ∼ N (µ, Σ) un vecteur gaussien (en vecteurs colonnes)
Y
avec ! !
µX ΣX ΣXY
µ= , Σ= ,
µY ΣY X ΣY
3. FILTRE DE KALMAN 21

La loi de X − E[X|Y ] est un loi gaussienne N (0, ΣX|Y ) avec

ΣX|Y = ΣX − ΣXY Σ−1

Y ΣY X .

Démonstration. Le caractère gaussien vient du fait que tout est combinaison linéaire d’un
vecteur gaussien. Ensuite, il suffit de calculer explicitement la covariance.

!
X
Résumé Dans le cadre des vecteurs gaussiens on peut écrire
Y

X = E(X|Y ) + X − E(X|Y )

avec E(X|Y ) et X − E(X|Y ) deux vecteurs gaussiens indépendants :

— E[X|Y ] = µX + ΣXY Σ−1
Y (Y − µY ), et en particulier, on retrouve que E(E[X|Y ]) =
µx .
— X − E(X|Y ) a pour loi N (0, ΣX|Y ) avec

ΣX|Y = ΣX − ΣXY Σ−1

Y ΣY X .

Comme corollaire, on peut donc obtenir des intervalles de confiance pour l’estima-
teur E[X|Y = y] de X.
Attention ! On dit parfois par abus de language, que la loi de X sachant Y est celle
d’un vecteur gaussien de moyenne E(X|Y ) et de matrice de covariance ΣX|Y .

3 Filtre de Kalman
On s’intéresse à un état xk ∈ Rn qui évolue selon un temps discret. On a donc une
suite (xk )k≥0 de variables aléatoires. On veut
— une estimation de l’état xk au cours du temps ;
— à partir d’observations incomplètes ou bruitées de xk ;
— connaissant un modèle théorique d’évolution pour xk .

Historiquement Ce problème s’est posé dans les années 1960 pour suivre la trajectoire
des modules Apollo. Il fallait estimer la position xk à partir des mesures terrestres yk
(soumis à une certaine incertitude), et connaissant un modèle physique xk+1 = F xk
(déplacement d’un point dans un champ de gravité) soumis encore une fois à une certaine
incertitude des forces environnantes.
22 CHAPITRE II. MODÈLE GAUSSIEN (2 SEMAINE)

Le modèle Précisons notre modèle théorique faisant intervenir un certain aléas. On a

un modèle linéaire pour l’évolution de (xk )k≥0

xk+1 = F xk + k

où F est une matrice n × n connue et (k )k≥0 une suite indépendantes de gaussiennes
centrée N (0, Q) de taille n.
Malheureusement on a ni accès aux (k )k≥0 ni aux (xk )k≥0 . Pour connaître xk , on a des
observations indirectes. On a accès à des mesures yk ∈ Rm qui sont des transformations
linéaires bruitées de xk :
yk = Hxk + ηk
où H est une matrice m × n et (ηk )k≥0 une suite indépendantes de gaussiennes centrée
N (0, R) de taille m.
Pour avoir un modèle théorique d’évolution, il est parfois utile d’intégrer dans l’état xk
des grandeurs auxiliaires en plus des grandeurs d’intérêt, comme dans les deux exemples
suivants, où on rajoute la vitesse en plus de la position.

Exemple 1 : Si on veut modéliser la position d’un véhicule avec deux coordonnées (x, y),
on utilise pour notre modèle une version discrétisée (avec un pas h = ∆t) de l’équation
position/vitesse v x = dtd x et v y = dtd y : c’est-à-dire xk+1 = xk + h · vkx et yk+1 = yk + h · vky .
Le modèle d’évolution des coordonnées positions/vitesses s’écrit alors
   
xk+1 xk
 y   y 
 k+1   k 
 x  = F ·  x  + k ,
 vk+1   vk 
y
vk+1 vky
avec  
1 0 h 0
 0 1 0 h 
F = .
 
 0 0 1 0 
0 0 0 1
Le GPS nous donne la position, donc on récupère les coordonnées bruitées
 
! x k
x̃k  y 
 k 
= H ·  x  + ηk ,
ỹk  vk 
vky
avec !
1 0 0 0
H= .
0 1 0 0
3. FILTRE DE KALMAN 23

Exemple 2 : Si on veut modéliser les oscillations x(t) d’un pendule soumis à l’équation

x00 (t) = −ω 2 x(t),

on va formuler
! l’évolution sous la forme!vectorielle y 0 (t) = Ay(t) pour le vecteur y(t) =
x(t) 0 1
0
et la matrice A = . On a donc y(t + h) = ehA y(t). Pour utiliser le
x (t) −ω 2 0
filtre de Kalman, on écrit notre modèle sous une version discrétisée (avec un pas h = ∆t) :

zk+1 = ehA zk + σk ,

!
σ2 0
pour des vecteurs gaussiens de R2 , avec σ > 0. La variance du bruit σk est .
0 σ2
On peut vérifier que !
1
cos(ωh) sin(ωh)
ehA = ω
−ω sin(ωh) cos(ωh)
(en dérivant par rapport à h par exemple, ou en diagonalisant A). On ne mesure que la
position, donc on récupère la première coordonnée de zn bruitée

wk = H · zk + δηk ,
 
avec H = 1 0 , et ηn des gaussiennes (la variance du bruit δηk est δ 2 ).

Idée générale du filtre de Kalman On va estimer l’état xk+1 à partir d’une estimation
de xk (calculée précédemment) et de l’observation yk+1 .
Cette estimation se fait en deux étapes :
— On appellera x̂k|k l’estimation de xk à partir des données jusqu’à l’instant k : c’est-
à-dire l’espérance conditionnelle de xk sachant (y0 , . . . , yk ).
— On va noter d’autre part x̂k+1|k l’estimation de xk+1 à partir des données jusqu’à
l’instant k : c’est-à-dire l’espérance conditionnelle de xk+1 sachant (y0 , . . . , yk ).
—
D’après la section précédente,

x̂k|k = E (xk |(y0 , . . . , yk ))

et


xk − x̂k|k ∼ N 0, Pk|k .
on dit alors que la loi de xk sachant (y0 , . . . , yk ) est une gaussienne N (x̂k|k , Pk|k ).
De façon similaire, on note

E(xk+1 |(y0 · · · , yk )) = x̂k+1|k

24 CHAPITRE II. MODÈLE GAUSSIEN (2 SEMAINE)

et
xk+1 − E(xk+1 |(y0 · · · , yk )) ∼ N (0, Pk+1|k )
Les espérances conditionnelles x̂k|k et x̂k+1|k et les matrices de covariance Pk|k et Pk+1|k
sont des quantités qu’on peut calculer au fur et à mesure grâce à la proposition suivante.

Proposition 3.1 (Filtre de Kalman). Récursivement, on a

— x̂k+1|k = F x̂k|k ;
— Pk+1|k = F · Pk|k · F T + Q ;
et d’autre part on notera Kk+1 = Pk+1|k H T (HPk+1|k H T + R)−1 , on a alors
— x̂k+1|k+1 = x̂k+1|k + Kk+1 (yk+1 − H x̂k+1|k ) ;
— Pk+1|k+1 = (I − Kk+1 H)Pk+1|k .

Comme prévu, on utilise donc les formules en 2 temps :

— Etape de prédiction pour passer de (x̂k|k , Pk|k ) à (x̂k+1|k , Pk+1|k ) ;
— Etape de correction pour passer de (x̂k+1|k , Pk+1|k ) à (x̂k+1|k+1 , Pk+1|k+1 ).
Remarquons que l’estimation initiale x̂0|0 et sa variance P0|0 sont nécessaires mais pas
toujours disponibles. Si l’état initial x0 est connu avec certitude, prendre x̂0|0 = x0 et
P0|0 = 0. Sinon, il faut utiliser une première mesure y0 de x0 .

Démonstration. Pour obtenir les formules remarquons que xk+1 = F xk + k . En condi-

tionnant par (y0 , . . . , yk ), on obtient

E(xk+1 |(y0 , . . . , yk )) = F E(xk |(y0 , . . . , yk )) + E(k |(y0 , . . . , yk ))

Comme k est indépendant de x0 , · · · xk et donc de y0 , · · · yk , E(k |(y0 , . . . , yk )) = E(k ) =

0. On en déduit donc que
x̂k+1|k = F x̂k|k
Pour ce qui est de la variance :

Pk+1|k = Var (xk+1 − E(xk+1 |(y0 · · · yk ))

= Var(xk+1 − x̂k+1|k )
= Var(F (xk − x̂k|k ) + k )

Or comme les deux variables alétaoires sont indépendants, leurs variances se somment :

Pk+1|k = Var(F (xk − x̂k|k )) + Var(k )

= F · Pk|k · F T + Q .

Maintenant, remarquons que

! ! !
xk+1 In 0
= · xk+1 + .
yk+1 H ηk+1
3. FILTRE DE KALMAN 25

En conditionnant par (y0 , . . . , yk ),

! ! !
xk+1 x̂k+1|k
E |(y0 , . . . , yk ) =
yk+1 H · x̂k+1|k

et
! ! !  ! !T !
xk+1 xk+1 In In 0 0
−E |(y0 , . . . , yk ) ∼ N 0, · Pk+1|k · + 
yk+1 yk+1 H H 0 R

Cette matrice de covariance se réécrit

!
Pk+1|k Pk+1|k · H T
H · Pk+1|k H · Pk+1|k H T + R

Pour en déduire la loi de xk+1 conditionnée par (y0 , . . . , yk+1 ), on travaille conditionnel-
lement à (y0 , · · · yk ). Il suffit alors de conditionner la première !
coordonnée par la seconde
xk+1
coordonnée dans cette loi conditionnelle gaussienne de sachant (y0 , · · · yk ). Dans
yk+1
ce cas là on peut montrer (admis ici) que le résultat de la section précédente sur le condi-
tionnement s’applique également et on obtient que l’espérance conditionnelle de xk+1 par
rapport à (y0 , . . . , yk+1 ) vaut

x̂k+1|k+1 = x̂k+1|k + Pk+1|k H T (H · Pk+1|k H T + R)−1 · (yk+1 − H x̂k+1|k )

et que l’écart xk+1 − x̂k+1|k+1 a pour variance

Pk+1|k+1 = Pk+1|k − Pk+1|k H T (H · Pk+1|k H T + R)−1 HPk+1|k ,

ce qui correspond exactement au formule annoncée dans la proposition.

On peut aussi réunir les deux étapes et calculer directement Pk+1|k en fonction de
Pk|k−1 . C’est ce qu’on appelle l’équation de Riccati, qui s’écrit alors

−1
Pk+1|k = Q + F Pk|k−1 F T − F Pk|k−1 H T · HPk|k−1 H T + R · HPk|k−1 F T .
26 CHAPITRE II. MODÈLE GAUSSIEN (2 SEMAINE)
Chapitre III

Processus de Poisson et Files d’attentes

(1 semaine)

Pré-requis : Lois de Poisson, Lois exponentielles.

Dans cette partie, on introduit le processus de Poisson qui est le modèle le plus simple
pour décrire les occurrences dans le temps d’un certain phénomène. Ce processus et ses
extensions peuvent-être utilisés dans de nombreux domaines : modéliser les passages suc-
cessifs d’un bus à un arrêt, les instants de pannes de composants électroniques, les arrivées
successives de clients dans un magasin, l’évolution de tailles de population...

1 Simulations de variables aléatoires de loi de Poisson

Pour rappel une variable aléatoire de Poisson de paramètre λ > 0 est une variable
aléatoire X à valeurs dans N = {0, 1, 2, · · · } telle que

λk
P(X = k) = exp(−λ), ∀k ∈ N.
k!
En particulier, E(X) = λ et P (X = 0) = exp(−λ).
Dans cette section, nous allons voir comment simuler de telles variables aléatoires.

1.1) Méthode récursive

La méthode de simulation des variables alétaoires discrète à partir d’une variable
uniforme est toujours possible :

Proposition 1.1. On pose s−1 = 0 et pour tout k ≥ 0,

k k
X X λi
sk = P (X = i) = exp(−λ)
i=0 i=0
i!

27
28 CHAPITRE III. POISSON ET FILES D’ATTENTES

Soit U est une variable aléatoire de loi uniforme U([0, 1]) et

∞
X
X= k1(sk−1 ≤U <sk ) .
k=0

Alors, X est une variable aléatoire de loi de Poisson de paramètre λ.

Toutefois cette méthode est difficile à implémenter telle quelle numériquement car elle
nécessite de connaitre toutes les sommes cumulées de probabilités sk . Un astuce consiste
en fait à calculer récursivement uniquement les probabilités cumulées nécessaire :

Simulation récursive d’une variable de Poisson

U =uniform()
s = exp(−λ)
p=s
x=0
while (U > s) :
x=x+1
p = p λx
s=s+p
return x;

-> cf TD7 : Exercice 2

1.2) A partir de variables exponentielles

Théorème 1.2. Soit (Xn )n≥1 une suite de variables aléatoires indépendantes de loi ex-
ponentielle de paramètre λ et
X n
Sn = Xi
i=1

avec par convention S0 = 0. Alors

X
N= 1{Sn ≤1} = max{n tels que Sn ≤ 1}
n≥1

suit une loi de Poisson de paramètre λ.

On déduit donc que l’on peut simuler une variable aléatoire de Poisson de la façon
suivante :
1. SIMULATIONS DE VARIABLES ALÉATOIRES DE LOI DE POISSON 29

Simulation exponentielle d’une variable de Poisson

x=0
S =variable aléatoire de paramètre λ
while (S ≤ 1) :
S = S+variable aléatoire de paramètre λ
x=x+1
sortir x;
Démonstration. Soit k ≥ 0, comme la suite (Sn ) est une suite croissante (comme somme
de variables aléatoires positives),
X
P(N = k) = P( 1{Sn ≤1} = k)
n≥1

= P(Sk ≤ 1 < Sk+1 )

Remarquons maintenant que

k+1
X k
X
Sk+1 = Xi = Xi + Xk+1 = Sk + Xk+1
i=1 i=1

et que par définition, les variables aléatoires Xk+1 et Sk sont indépendantes. Aussi si l’on
note fSk et fXk+1 leurs densités respectives, on a
Z Z
P (Sk ≤ 1 < Sk+1 ) = 1s≤1<x+s fSk (s)fXk+1 (x)dsdx.
Z Z
= 1s≤1 1x>1−s fSk (s)fXk+1 (x)dsdx.
Z 1 Z ∞ 
= dsfSk (s) dxfXk+1 (x)
0 1−s

Rappelons maintenant que Sk est la somme de k variables aléatoires exponentielles indé-

pendantes, aussi Sk suit une loi gamma de paramètres (λ, k) et la densité
exp(−λs)
fSk (s) = λk sk−1 .
(k − 1)!
On en déduit que
Z 1 Z ∞ Z 1 Z ∞
k k−1 exp(−λs)
dsfSk (s) dxfXk+1 (x) = dsλ s λ exp(−λ(x))dx
0 1−s 0 (k − 1)! 1−s
Z 1
exp(−λs)
= dsλk sk−1 exp(−(1 − s)λ)
0 (k − 1)!
exp(−λ)
= λk
k!
Ce qui conclut la preuve.
30 CHAPITRE III. POISSON ET FILES D’ATTENTES

2 Processus de Poisson
Nous allons maintenant introduire le processus de Poisson. C’est un premier exemple
de processus en temps continu.

Définition 2.1. Soit (Xn ) une suite de variables aléatoires i.i.d. de lois exponentielles
de paramètre λ, on définit comme précédemment Sn = nk=1 Xk et S0 = 0. Pour chaque
P

t ≥ 0 on définit un processus de Poisson de paramètre λ par

X
Nt = 1{Sn ≤t} = max{n ≥ 0 tels que Sn ≤ t}.
n≥1

La définition ci-dessous correspond en fait à la construction du processus et non à sa

vraie définition (mais c’est amplement suffisant).
Le processus de Poisson est un processus de comptage, c’est à dire qu’il compte le
nombre d’évènements qui se sont produits dans une fenêtre de temps [0, t].

Quelques propriétés

Proposition 2.2. Soit Nt un processus de Poisson de paramètre λ > 0, alors

i) La fonction t 7→ Nt est croissante, elle part de N0 = 0 et prend ses valeurs dans
N.
ii) Quelque soit 0 ≤ t la variable aléatoire Nt est une variable aléatoire de Poisson de
paramètre λt.
Plus généralement pour tout 0 ≤ a ≤ b, Nb − Na est une variable aléatoire de
Poisson de paramètre λ(b − a) qui représente le nombre de points dans l’intervalle
]a, b].
iiii) Conditionnellement à l’évènement {Nt = n}, la loi des instants de saut (S1 , · · · Sn )
est la loi d’un échantillon de n variables aléatoires uniformes sur [0, t] réordonnées
par ordre croissant. C’est à dire qu’elle a la densité
n!
10≤t1 <t2 <tn ≤t
tn
Démonstration. i) Pour tout n fixé 1Sn ≤t est une fonction croissante, en effet elle vaut 0
si Sn > t (donc pour t petit), puis 1 lorsque t dépasse Sn .
ii) La preuve est très similaire à celle ci-dessus (en exercice).

P(Nt = n) = P(Tn ≤ t < Tn+1 )

(λt)n
= exp(−λt)
n!
2. PROCESSUS DE POISSON 31

iii) Commençons par calculer la densité de (S1 , · · · Sn ).

E(f (S1 , · · · , Sn )) = E(f (X1 , X1 + X2 , · · · , X1 + .... + Xn ))

Z
= f (x1 , x1 + x2 , · · · , x1 + .... + xn ))λn exp(−λx1 )... exp(−λ(xn ))dx1 · · · dxn
Z
= f (x1 , t2 , · · · tn )λn exp(−λtn )10≤x1 <t2 <...<tn dt1 · · · dtn

Maintenant pour tout ensemble mesurable Γ ⊂ Rn , on a

P((S1 , · · · Sn ) ∈ Γ, Nt = n)
P((S1 , · · · Sn ) ∈ Γ|Nt = n) =
P(Nt = n)
Z
n!
= nn 1sn ≤t<sn+1 1(s1 ,...,sn )∈Γ λn+1 exp(−λsn+1 )10≤s1 <s2 <...<sn <sn+1 ds1 · · · dsn+1
λ t exp(λt)
Z Z ∞
n!λ
= n 10≤s1 <s2 <...<sn <t exp(−λsn+1 )dsn+1
t Γ t
Z
n!
= n 10≤s1 <s2 <...<sn <t ds1 ...dsn
t Γ

Finalement il reste à montrer que si on considère un tirage de (U1 , ...Un ) variables aléatoires
uniformes indépendantes sur [0, t], alors le n-uplet trié (U(1) , ...U(n) ) a la même densité.
Ceci provient simplement du fait qu’il existe n! ordre pour ces tirages et que tous ont la
même probabilité de ce produire, ce qui donne facilement le résultat.

On déduit de cette proposition une seconde methode pour simuler un processus de

Poisson sur un intervalle de temps [0, T ] :
— On tire le nombre de point NT suivant une loi de Poisson de paramètre λT
— Sachant la valeur de NT , on tire variables aléatoires uniformes sur [0, T ] et on les
ordonne.
Le processus de Poisson se caractérise aussi par de bonnes propriétés d’indépendances

Proposition 2.3. Considérons deux intervalles de temps disjoints ]T1 , T2 ] et ]T3 , T4 ] alors
les variables aléatoires NT2 − NT1 et NT4 − NT3 sont des variables aléatoires indépendantes
de lois de Poisson de paramètres respectifs (T2 − T1 )λ et (T4 − T3 )λ

Comportement en temps long

Proposition 2.4. Soit (Nt )t≥0 un processus de Poisson de paramètre λ > 0, alors

Nt
−→ λ p.s..
t t→∞
32 CHAPITRE III. POISSON ET FILES D’ATTENTES

Démonstration. Commençons par écrire pour un entier n

n
X
Nn = Nk − Nk−1
k=1

Ainsi, comme les variables aléatoires Nk − Nk−1 sont indépendantes et de même lois
P oiss(λ), on peut appliquer la loi forte des grands nombre pour en déduire que
Nn
−→ λ p.s..
n n→∞
Ensuite, on utilise la croissance du processus Nt pour en déduire que pour tout t ∈
[n, n + 1],
Nn n Nt Nn+1 n + 1
≤ ≤
n+1n t n n+1
En faisant tendre n → ∞, on obtient la conclusion.

Il existe de nombreuses généralisation des processus de Poisson, que ce soit dans des
espaces spatiaux temporels, ou avec des fonctions d’intensités qui dépendent du temps.

3 Files d’attentes
Dans la suite, nous allons voir une utilisation des processus de Poisson à la modélisation
des files d’attentes. L’étude des files d’attente est un domaine actif très utile notamment
pour la gestion des serveurs, des requêtes et dans de nombreux domaines d’ingénierie.
Une file d’attente est caractérisée par trois paramètres
— la loi décrivant les arrivées successives de clients dans la file. On utilise souvent
un processus de Poisson, c’est à dire que l’on suppose qu’entre deux arrivées de
clients, il se passe un temps exponentiel de paramètre λ.
— une loi décrivant le temps nécessaire pour servir un client. Ce sont des variables
aléatoires que l’on considère en général comme indépendantes et de même loi. Là
encore, c’est la loi exponentielle qui donne le traitement mathématique le plus
simple.
— un paramètre qui décrit le nombre de serveurs disponibles : ce peut être un entier
N > 0 qui correspond au nombre de guichet ou bien l’infini. Lorsque le nombre de
serveur est infini chaque client qui arrive est immédiatement pris en charge.

3.1) File d’attente M/M/1

L’exemple le plus simple de file d’attente s’appelle la file M/M/1 dans laquelle, il y
a un unique serveur, les clients arrivent selon un processus de Poisson de paramètre λ et
3. FILES D’ATTENTES 33

sont servis en un temps exponentiel de paramètre β. Le ’M’ représente alors le fait que la
loi exponentielle est sans mémoire et donc la file d’attente est un processus de Markov à
temps continu.
On souhaite étudier le nombre Wt de personnes dans la file d’attente. C’est la différence
entre le nombre d’arrivées At et le nombre de départ Dt . Ainsi, lorsque Wt = 0, la file
comme le guichet sont libres, lorsqueWt = 1 ceci signifie qu’une personne est au guichet,
et personne en attente, et lorsque Wt ≥ 2, le guichet est occupé et Wt − 1 personnes sont
en attentes.

Rappelons des propriétés importantes des variables exponentielles :

Proposition 3.1. Soit A et B deux variables exponentielles de paramètres λ et β, alors
— V = min(A, B) est une variable exponentielle de paramètre λ + β
— D’autre part, la variable aléatoire W = 1V =A suit une loi de Bernoulli de paramètre
λ
λ+β
.
— V et W sont indépendantes.
Ceci nous permet en particulier de simuler la file d’attente de la façon suivante. Partons
d’une nombre n de clients dans la file d’attente W0 = n.
— Si n = 0 le prochain évènement est l’arrivée d’un client qui se produit après un
temps exponentiel de paramètre λ.
— Si n > 0, le prochain évènement se produit selon une loi exponentielle de paramètre
λ
λ + β. Avec probabilité λ+β on rajoute un nouveau client, et sinon on en enlève 1.
Le mieux pour simuler une file d’attente est de construire deux vecteurs. Le premier vec-
teur contient les temps successifs Tn des différents évènements, le second vecteur contient
le nombre de clients WTn dans la file à ces différents instants.

Chaine de Markov incluse On peut choisir de n’étudier que les variations du nombre
de clients dans la file sans prendre en compte l’aspect temporel. On définit alors un
processus à temps discret Zn = WTn où Tn est le temps aléatoire du nme évènement.
On peut montrer que Zn est une chaine de Markov sur N dont les transitions sont
décrites par

Q(0, 1) = 1
λ
Q(i, i + 1) = = 1 − Q(i, i − 1), i≥1
λ+β
On a alors des résultats sur le comportement asymptotique de la chaine de Markov
(Zn )n≥0 .
Pour un processus irréductible, comme pour une chaîne irréductible apériodique, trois
cas sont possibles en fonction du retour dans un état donné.
34 CHAPITRE III. POISSON ET FILES D’ATTENTES

— Si le processus, partant d’un état donné, a une probabilité non nulle de ne jamais y
retourner, il est dit transient. En pratique cela signifie que le nombre que l’on étudie
(taille de population ou nombre de clients) tend vers l’infini. Nous interprétons cette
situation comme une saturation du système.
— Si le processus, partant d’un état donné, y revient nécessairement au bout d’un
temps fini en moyenne, il est dit récurrent positif. Dans ce cas, il converge en loi
vers sa mesure stationnaire, interprétée comme une situation d’équilibre du système
étudié.
— Si le processus, partant d’un état donné, y revient avec probabilité 1 mais que
l’espérance du temps de retour est infinie, il est dit récurrent nul.
On peut alors montrer le résultat suivant :

Proposition 3.2. La chaine Zn est

— transiente si λ > β
— récurrente positive si λ < β
— récurrente nulle si λ = β

Démonstration. Voir en TD

Comportement limite de la chaine en temps continu On s’interesse à ce qui se

passe lorsque λ < β pour le processus en temps continu. On admettra ici les résultats
suivants
(loi)
Proposition 3.3. Lorsque λ < β, la file d’attente Wt =⇒ π où π est la mesure station-
naire donnée par
   n
λ λ
πn = 1 − .
β β

En particulier, on peut déduire de ce résultat plusieurs comportement de la file d’at-

tente en régime stationnaire :
— le nombre moyen de clients dans la file d’attente est donné par

λ
E(π) = .
β−λ

— la loi du temps passé dans la file d’attente par un client est donnée en régime
stationnaire par E(β − λ).
Pour ce dernier résultats, remarquons que quand la chaine est vide, le temps d’attente
d’un client est E(β).
Lorsqu’il y a déjà n clients devant lui, son temps d’attente est donc la somme des n + 1
3. FILES D’ATTENTES 35

services des n clients précédents et du sien. Comme chaque service est i.i.d de loi E(β), la
loi du temps de service a pour densité

β n+1 tn
exp(−βt)
n!
Si le nombre de client suit la loi π on en déduit donc que la densité du temps de service
est
X β n+1 tn X  λ n β n+1 tn
fT (t) = πn exp(−βt) = exp(−βt)
n≥0
n! n≥0
β n!
 
λ
=β 1− exp(−βt) exp(λt)
β
= (β − λ) exp(−(β − λ)t).

3.2) D’autres modèles

Il existe des modèles plus généraux qui prennent en compte différents serveurs, un
nombre maximum de clients dans la file, des lois de services non exponentielles.... Nous
verrons en TP comment simuler ces files d’attentes.
36 CHAPITRE III. POISSON ET FILES D’ATTENTES
Chapitre IV

Branchement (1 semaine)

Pré-requis : Martingales, espérance conditionnelle

1 Fonction génératrice
Soit X une variable aléatoire discrète sur N, de loi déterminée par les probabilités

pX (k) = P[X ∈ k], ∀k ∈ N.

Définition 1.1. La fonction génératrice de X est la fonction GX : [−1, 1] → [−1, 1]

définie par
X
GX (s) = E[sX ] = pX (k) · sk .
k≥0

Exemples :
— si X suit la loi de Poisson P(λ) de paramètre λ > 0. Alors
X λk
GX (k) = e−λ · sk = esλ e−λ = e(s−1)λ .
k≥0
k!

— si X suit la loi Binomiale B(n, p) de paramètre n ≥ 1, 0 < p < 1. Alors

n  
X n  n
GX (k) = pk (1 − p)n−k · sk = (1 − p) + ps .
k=0
p

— si X suit la loi géométrique G(p) de paramètre 0 < p < 1. Alors

X sp
GX (k) = (1 − p)k−1 p · sk = .
k≥1
1 − (1 − p)s

Proposition 1.2. La fonction génératrice GX satisfait les propriétés suivantes :

37
38 CHAPITRE IV. BRANCHEMENT (1 SEMAINE)

1. elle est analytique sur ]−1, 1[ (c’est-à-dire développable en série entière au voisinage
de chacun des points de ] − 1, 1[).
2. elle est croissante sur [0, 1].
3. elle est convexe (et même strictement convexe si P[X ≥ 2] > 0) sur [0, 1].
1
4. elle caractérise la loi de X car pX (n) = n!
G(n) (0).
5. si E[X] < +∞, alors GX est dérivable (à gauche) en 1 et G0X (1) = E[X].
6. GX (0) = pX (0) = P[X = 0].

Démonstration. 1. GX (1) = 1 < +∞ donc la série entière a un rayon de convergence

supérieure ou égale à 1.
2. pour s ∈ [0, 1[, on a donc
X
G0 (s) = kpX (k) · sk−1 ≥ 0
k≥1

3. pour s ∈ [0, 1[, on a donc

X
G00 (s) = k(k − 1)pX (k) · sk−2 ≥ 0
k≥2

et même G00 (s) > 0 dès que s ∈]0, 1[ et qu’un des pX (k) est strictement positif pour
k ≥ 2.
4. Plus généralement, sur ] − 1, 1[,
X
G(n) (s) = k(k − 1) · · · (k − n + 1)pX (k) · sk−n
k≥n

et donc seul le terme k = n compte pour calculer G(n) (0) = n! · pX (n).

5. on a G0 (s) = k≥1 kpX (k) · sk−1 pour s ∈] − 1, 1[ et cette somme converge vers
P

E[X] lorsque s converge vers 1 (par convergence monotone par exemple). On déduit
alors du théorème des accroissements finis que cette limite est la dérivée à gauche
G0X (1).
6. Par définition, GX (0) = E[sX ] = E[1X=0 ] = pX (0).

Les fonctions génératrice sont particulièrement adapté pour traiter la somme de va-
riables aléatoires discrètes indépendantes.

Proposition 1.3. Si X et Y sont des variables aléatoires indépendantes à valeurs dans

N, on a GX+Y = GX GY .
2. PROCESSUS DE GALTON-WATSON 39

Démonstration. On écrit E[sX+Y ] = E[sX eY ] = E[sX ]E[sY ].

Plus généralement, il est possible de traiter une somme d’un nombre aléatoire de
variables aléatoires.

Théorème 1.4. Soit N une variable discrète sur N et (Xn )n≥1 des variables aléatoires
indépendantes à valeurs dans N, de même loi et donc de même fonction génératrice GX ,
et indépendantes de N . Alors, considérant la variable aléatoire
N
X
SN = Xk ,
k=1

on a GSN = GN ◦ GX , au sens où, pour s ∈ [−1, 1],

GSN (s) = GN (GX (s)).

P
Démonstration. L’idée est de décomposer selon les valeurs de N . On a 1 = n 1N =n donc
∞ ∞ h PN i X∞ h Pn i
E 1N =n s( k=1 Xk ) = E 1N =n s( k=1 Xk )
X X
E[sSN ] = E[1N =n sSN ] =
n=0 n=0 n=0

mais comme
h Pn i h Pn i n
E 1N =n s( k=1 Xk ) = P[N = n]E s( k=1 Xk ) = P[N = n]
Y
E sXk = P[N = n]GX (s)n ,
 
k=1

on obtient
∞
X
SN
E[s ]= P[N = n]GX (s)n = GN (GX (s)).
n=0

2 Processus de Galton-Watson
On considère une population d’individus (ce peut être n’importe quoi, des particules,
des gènes, des animaux etc) dont la taille évolue au cours d’un temps discret. Pour encoder
la taille, on va donc considérer une suite de variables aléatoires (Xn )n≥0 discrètes sur N.
Ces variables ne sont pas indépendantes, puisque le nombre d’invidus à un certain instant
n + 1 dépend du nombre d’individus à l’instant d’avant n avec une règle simple : chaque
individu va engendrer k individus avec la probabilité pk , indépendamment des autres
individus. On a donc le modèle suivant.
40 CHAPITRE IV. BRANCHEMENT (1 SEMAINE)

Définition 2.1. Soit (Yn,k )n,k∈N une suite de variables aléatoires indépendantes de même
loi. On définit le processus de Galton-Watson comme la suite de variables aléatoires
(Xn )n≥0 de la façon suivante : X0 = 1 et pour n ≥ 0,
Xn
X
Xn+1 = Yn,k .
k=1

Dans notre modèle, il y a un seul individu à la génération 0. Chaque individu k de la

n-ième génération a un nombre Yn,k de descendants (pour 1 ≥ k ≥ Xn ). Ce nombre de
descendant suit toujours la même loi et est indépendant des autres. On va appeler G la
fonction génératrice commune des Yn,k , et on va supposer que la moyenne m = E[Yn,k ] est
finie, de telle sorte que
G0 (1) = m,
et on va supposer que la probabilité de ne pas produire de descendant est strictement
positive :
G(0) > 0.
Le théorème de la section précédente nous dit que

GXn+1 = GXn ◦ G

puisque Xn est fonction des (Ym,k )m<n,k et est donc indépendant des (Yn,k )k . De cette
formule, on peut déduire que :

Proposition 2.2. La taille moyenne du processus est E[Xn ] = mn .

Démonstration. Il suffit de remarquer que

E[Xn+1 ] = G0Xn+1 (1) = G0 (1)G0Xn (G(1)) = m · G0Xn (1) = mE[Xn ]

et que par conséquent, (E[Xn ])n est une suite géométrique de raison m.

Intéressons-nous maintenant à la probabilité d’extinction de la population. Remar-

quons que l’événement qui correspond au fait que la population s’éteigne au bout d’un
moment s’écrit ∪n {Xn = 0}. C’est une union croissante d’événements et la probabilité
q correspondant à l’extinction de la population est donc la limite croissante de la suite
xn := P[Xn = 0], c’est-à-dire
q = lim xn .
n→∞

On a xn = P[Xn = 0] = GXn (0), et si on veut utiliser la formule GXn+1 = GXn ◦ G pour

calculer xn , il faut la reformuler.

Proposition 2.3. Pour n ≥ 1, GXn = |G ◦ ·{z

· · ◦ G}.
n fois
2. PROCESSUS DE GALTON-WATSON 41

Démonstration. On a d’abord GX1 = GY0,1 = G, puis par récurrence, si

GXn = |G ◦ ·{z
· · ◦ G}
n fois

alors
GXn+1 = GXn ◦ G = G
| ◦ ·{z
· · ◦ G} ◦G = G
| ◦ ·{z
· · ◦ G} .
n fois n+1 fois

Une conséquence directe est que

· · ◦ G}(0) = G ◦ G
xn+1 = GXn+1 = |G ◦ ·{z | ◦ ·{z
· · ◦ G}(0) = G(GXn (0)) = G(xn ).
n+1 fois n fois

La limite de xn est donc nécessairement un point fixe de G sur [0, 1] : G(q) = q. De la

formule xn+1 = G(xn ), on peut aussi en déduire que q est en fait le plus petit point fixe
sur [0, 1] car si q̃ est un autre point fixe de G sur [0, 1], on a forcément xn ≤ q̃ pour tout
n ≥ 1 : ça vient du fait que x1 = G(0) ≤ G(q̃) = q̃, et par récurrence xn ≤ q̃ entraîne que
xn+1 = G(xn ) ≤ G(q̃) = q̃.
Chercher la probabilité d’extinction q de (Xn )n≥0 , c’est donc chercher le plus petit
point fixe de G sur [0, 1]. La proposition suivante distingue deux cas.

Proposition 2.4. Si m = G0 (1) ≤ 1, alors q = 1 : la population s’éteint presque sûrement.

Si m = G0 (1) > 1, alors q < 1 : la population a une probabilité non-nulle de survivre.

Démonstration. La série génératrice G s’écrit

X
G(s) = pk · sk
k≥0

P
avec G(1) = k≥0 pk = 1. On va étudier la fonction F (s) = G(s) − s, qui va de F (0) =
G(0) > 0 à F (1) = 1 − 1 = 0. En effet, les points fixes de G sont les réels où F s’annule.
Sa dérivée est F 0 (s) = G0 (s) − 1 sur [0, 1], et on remarque que F 00 (s) = G00 (s) ≥ 0 donc
cette dérivée F 0 (s) est croissante (strictement si k≥2 pk > 0).
P

Supposons d’abord que m = G0 (1) ≤ 1.

P
Ou bien k≥2 pk = 0, et F (s) est alors une fonction linéaire qui décroît de F (0) > 0
à F (1) = 0 et qui ne s’annule donc qu’en s = 1 : le seul point fixe de G est q = 1.
Ou bien k≥2 pk > 0 et F 0 est strictement croissante. Comme F 0 (1) = G0 (1) − 1 ≤ 0,
P

on a F 0 (s) < 0 pour s ∈ [0, 1[ ce qui signifie que F (s) est strictement décroissante de
F (0) > 0 à F (1) = 0. Elle ne s’annule donc qu’en s = 1 : le seul point fixe de G est q = 1.
Supposons maintenant que m = G0 (1) > 1.
42 CHAPITRE IV. BRANCHEMENT (1 SEMAINE)

Comme F 0 (1) > 0, la fonction F 0 est strictement positive dans un voisinage de 1 et

donc F est strictement croissante dans un voisinage de 1. Il existe donc s dans un voisinage
autour de 1 tel que F (s) < F (1) = 0. Mais comme F (0) > 0, le théorème des valeurs
intermédiaire nous dit qu’il existe alors un réel q̃ ∈]0, s[ tel que F (q̃) = 0 : il existe un
point fixe de G plus petit que 1, donc q < 1.

3 Cas sous-critique
On a vu que si m ≤ 1, la population s’éteint presque sûrement. Pour autant, il y a une
différence entre le cas m < 1 et le cas m = 1. On peut par exemple regarder l’espérance
de la taille totale de la population
"∞ # ∞ ∞
X X X
E Zn = E [Zn ] = mn
n=0 n=0 n=0

qui est finie si m < 1 et qui est infinie si m = +∞. Le théorème suivant permet de
quantifier la convergence de P[Xn > 0] vers 0 lorsque n grandit.

Théorème 3.1. Si m = E[Yn,k ] < 1, alors

P[Xn > 0] ≤ mn .

Démonstration. Puisque E[Xn ] = E[Xn |Xn > 0]P[Xn > 0] + 0 × P[Xn = 0], On a

E[Xn ]
P[Xn > 0] = ≤ E[Xn ] = mn .
E[Xn |Xn > 0]

4 Cas critique
Si m = 1, l’espérance E[Xn ] reste toujours égale à 1. La taille de Xn lorsque la
population n’est pas éteinte doit donc être très grande pour maintenir une espérance
constante. En fait, on peut montrer que supn Xn n’est pas intégrable, car sinon, on aurait
par convergence dominée l’égalité absurde suivante

1 = lim E[Xn ] = E[lim Xn ] = 0.

n n

Autrement dit, la martingale Xn converge vers 0 p.s. mais n’est pas convergente au sens
L2 . Le théorème suivant permet de mieux comprendre la taille de Xn lorsque la population
n’est pas éteinte dans le cas critique.
4. CAS CRITIQUE 43

Théorème 4.1. Si m = E[Yn,k ] = 1 et σ 2 = Var[Yn,k

2
] < +∞, alors
1. limn nP[Xn > 0] = 2/σ 2 ;
2. limn n1 E[Xn |Xn > 0] = σ 2 /2.
Démonstration. En dérivant deux fois G, on obtient G00 (s) = k≥2 k(k − 1)p(k)sk−2 qui
P

converge vers E[Yn,k (Yn,k − 1)] = σ 2 lorsque s tend vers 1 (convergence monotone). Par
conséquent, G est deux fois dérivable à gauche en 1, et on a la formule de Taylor avec
reste intégral en 1 à l’ordre 2
1
G(s) = s + (s − 1)2 σ 2 + (s − 1)2 R(s)
2
ave R(s) borné sur [0, 1] et qui tend vers 0 quand s → 1. Autrement dit, on a

G(s) − s = (1 − s)2 h(s)

avec h(s) = σ 2 /2 + R(s). On en déduit que

1 1 G(s) − s
− =
1 − G(s) 1 − s (1 − G(s))(1 − s)
h(s)(1 − s)
=
1 − s − (1 − s)2 h(s)
h(s)
=
1 − (1 − s)h(s)
(1 − s)h(s)
= h(s) +
1 − (1 − s)h(s)
= σ 2 /2 + k(s)
h(s)
avec k(s) = R(s) + (1 − s) 1−(1−s)h(s) qui est donc bornée et qui tend vers 0 quand s → 1.
Par sommage téléscopique, on peut écrire
  n−1  
1 1 1X 1 1
−1 = −
n 1 − GXn (0) n k=0 1 − GXn+1 (0) 1 − GXn (0)
comme la moyenne de Césaro de la suite
 
1 1
− = σ 2 /2 + k(GXn (0))
1 − G(GXn (0)) 1 − GXn (0)
qui tend donc vers la limite σ 2 /2 puisque GXn (0) = P[Xn = 0] tend vers 1. On a donc
   −1
1 1 1
nP[Xn > 0] = n(1 − GXn (0)) = −1 +
n 1 − GXn (0) n
qui tend vers 2/σ 2 puis par conséquent
1 E[Xn ]
E[Xn |Xn > 0] =
n nP[Xn > 0]
qui tend vers σ 2 /2.
44 CHAPITRE IV. BRANCHEMENT (1 SEMAINE)

5 Cas sur-critique
Rappelons que si m > 1, la probabilité d’extinction q est l’unique point fixe < 1 de la
fonction génératrice associée à la loi de Yn,k . Il y a donc une probabilité non-triviale que
la population survive.
Pour étudier plus finement (Xn )n≥0 , on introduit la suite (Xn /mn )n≥0 . Son utilité vient
du fait que c’est une martingale pour la filtration naturelle Fn = σ(X0 , . . . , Xn ) :
 
Xn+1 E [ Xn+1 | Xn ]
E n+1
Fn =
m mn+1
P`
car Xn est une chaîne de Markov ; or E [Xn+1 | Xn = `] = k=1 E[Yn,k ] = `m donc
E [Xn+1 | Xn ] = Xn m et
 
Xn+1 E [Xn+1 | Xn ] mXn Xn
E n+1
Fn = n+1
= n+1 = n .
m m m m

Comme c’est une martingale positive, elle converge p.s. vers une variable aléatoire réelle
positive. En fait, contrairement au cas où m = 1, on a même ici la convergence L2 .

Théorème 5.1. Si m = E[Yn,k ] > 1 et σ 2 = Var[Yn,k

2
] < +∞, alors
1. la convergence de (Xn /mn )n≥0 vers L a lieu en moyenne quadratique ;
2. L est de moyenne 1 et de variance σ 2 /(m2 − m) ;
3. P[L = 0] = q.

Le dernier point signifie en fait que l’événement {L > 0} est égal presque sûrement à
l’événement de la survie de la population.

Démonstration. On vérifie que (Xn /mn )n≥0 est une martingale bornée dans L2 en calcu-
lant la variance
n−1
n 2
 X
E[(Xk+1 /mk+1 − Xk /mk )2 ].

E (Xn /m − 1) =
k=0

On calcule
Xn
n+1 n 1 X
Xn+1 /m − Xn /m = (Yn,k − m)
mn+1 k=1

ce qui donne
 !2 
`
1 X 1
E[(Xk+1 /mk+1 − Xk /mk )2 |Xn = `] = 2n+2 E  (Yn,k − m)  = 2n+2 `σ 2
m k=1
m
5. CAS SUR-CRITIQUE 45

par indépendance des Yn,k − m, donc

1
E[(Xk+1 /mk+1 − Xk /mk )2 ] = mσ 2
m2n+2
puis
n−1
 n 2
 X 1 2 1 σ2
E (Xn /m − 1) = σ →n→∞ .
k=0
mk+2 m2 1 − 1/m
Le deuxième point est une conséquence de la moyenne constante E[Xn /mn ] = 1 et de la
limite σ 2 /(m2 − m) de la variance de Xn /mn .
Pour le dernier point, on va conditionner par rapport à X1 :
∞
X ∞
X
P[L = 0] = P[L = 0|X1 = k]P[X1 = k] = P[L = 0]k P[X1 = k] = G(P[L = 0]).
k=0 k=0

On a P[L = 0] comme point fixe de G, donc P[L = 0] = q ou 1. Le deuxième cas étant

impossible puisque E[L] = 1, on obtient P[L = 0] = q.
46 CHAPITRE IV. BRANCHEMENT (1 SEMAINE)
Chapitre V

Méthode de Monte-Carlo par chaînes

de Markov (1,5 semaines)

Pré-requis : Chaines de Markov et théorèmes de convergence

Les méthodes de Monte-Carlo par chaînes de Markov proposent des estimateurs en uti-
lisant des propriétés de convergence de chaînes de Markov vu dans le cours de probabilités,
comme par exemple la loi des grands nombres.

1 Algorithme de Metropolis-Hastings
Le but de l’algorithme de Metropolis-Hastings est de simuler une mesure de probabilité
π sur un espace d’état E en simulant une chaîne de Markov dont la mesure invariante
est π. Sous de bonnes hypothèses, cette chaîne de Markov convergera donc en loi vers π
(ergodicité). Cet algorithme est en particulier utilisé dans le cas, où la mesure π n’est pas
connue exactement, mais à constante près.
Le cadre général est donc le suivant : la mesure de probabilité π s’écrit

f (x)λ(dx)
π(dx) = R ,
E
f (x)λ(dx)

où
— f est une fonction mesurable de E 7→ R+
— λ est une mesure positive de référence (la mesure de Lebesgue par exemple)
R
La constante de renormalisation E f (x)λ(dx) n’est pas connue. On écrit en général que
la mesure π est proportionnelle à f (x)λ(dx) :

π ∼ f (x)λ(dx).

47
48 CHAPITRE V. MCMC

1.1) Espace d’état fini

Soit π une mesure de probabilité sur un espace d’état fini E.
Trouver une chaîne de Markov dont la mesure invariante est π n’est pas difficile : par
exemple, on peut prendre X1 , . . . , Xn indépendants de loi π. Cette solution triviale ne
nous intéresse pas. Notre but est de trouver une chaîne de Markov (Xn )n≥0 qu’on sait
simuler même dans les cas où on ne sait pas simuler π.
Partons du principe qu’on sache simuler simplement une chaîne de Markov auxiliaire
(dont la mesure invariante n’est pas forcément π) d’après une matrice de transition Q :
E × E → [0, 1]. L’algorithme de Metropolis-Hastings modifie Q pour obtenir une chaîne
de Markov de mesure invariante π.

Définition 1.1. La chaîne de Markov (Xn )n≥0 définie par Metropolis-Hastings à partir
d’une mesure π sur E et d’une matrice de transition Q est la chaîne de Markov homogène
de matrice de transition
(  
π(y)Q(y,x)
Q(x, y) min 1, π(x)Q(x,y) si x 6= y
P (x, y) = P
1 − z6=x P (x, z) si x = y

La probabilité P (x, y) d’aller de x à y (pour y 6= x) est donc la probabilité

 d’aller de 
π(y)Q(y,x)
x à y en suivant la matrice de transition Q multipliée par un facteur min 1, π(x)Q(x,y)
que l’on appelle probabilité d’acceptation rejet.
On peut donc simuler un saut partant de x comme ceci : on choisit un candidat y
sur lequel sauter suivant la loi donnée par Q(x, ·), et on saute effectivement sur y avec
π(y)Q(y,x)
probabilité min 1, π(x)Q(x,y) , sinon on reste sur x.
Remarquons que la probabilité de garder ou de rejeter le saut fait intervenir le quotient
π(x)/π(y). En particulier, il suffit de connaître π sur E à constante près, ce qui s’avère
très important dans certains cas.
La procédure peut être retranscrite comme ceci.

Metropolis-Hastings
X0 = simulation selon une loi initiale quelconque ;
pour i entre 1 et n,
y = tirage alatoire selon la loi Q(Xi−1 , ·)
π(y)Q(y,Xi−1 )
si unif orm() < π(X i−1 )Q(Xi−1 ,y)

Xi = y
sinon faire
Xi = Xi−1
sortir (X0 , . . . , Xn );
1. ALGORITHME DE METROPOLIS-HASTINGS 49

Vérification On veut vérifier que la chaine Xn construite ci-dessus admet bien P

comme matrice de transition.
Précisons un peu les notations. A chaque étape i, on notera Yi la variable aléatoire
tirée selon la loi Q(Xi−1 , ·) et Ui la variable aléatoire uniforme utilisée dans le choix. En
particulier, Ui est indépendante et Yi et donc de Xi−1 .
Pour tout couple (x, y) ∈ E avec x 6= y
 
π(Yn+1 )Q(Yn+1 , Xn )
P(Xn+1 = y|Xn = x) = P(Yn+1 = y et Un+1 ≤ |Xn = x)
π(Xn )Q(Xn , Yn+1 )
 
π(y)Q(yx)
= P(Yn+1 = y et Un+1 ≤ |Xn = x)
π(x)Q(x, y)
  
π(y)Q(yx)
= min 1, Q(x, y)
π(x)Q(x, y)

où l’on a utilisé le fait que Un+1 est indépendant de Xn et Yn+1 .

Justification de l’algorithme Montrons que la mesure π est bien la mesure invariante

associée à la chaine de Markov de matrice P .

Proposition 1.2. On suppose que π(x) 6= 0 pour tout x ∈ E et que Q(x, y) > 0 ⇔
Q(y, x) > 0.
i) π est la mesure de probabilité invariante pour la chaîne (Xn )n≥0 de matrice de
transition P .
ii) Si de plus Q est irréductible, alors P l’est aussi et on a donc les convergences
usuelles vers π quand n → ∞.

Démonstration. Preuve de i) : On commence par montrer que la mesure π est réversible

pour P c’est à dire que pour tout couple (x, y)µ(x)P (x, y) = µ(y)P (y, x).
Commençons par vérifier qu’une mesure réversible est invariante. Remarquons que
X
µP (y) = µ(x)P (x, y)
x∈E
X
= µ(y)P (y, x)
x∈E
X
= µ(y) P (y, x) = µ(y).
x∈E

Vérifions maintenant cette propriété. La propriété est toujours vraie si x = y.

50 CHAPITRE V. MCMC

Pour x 6= y, on a π(x)P (x, y) = 0 = π(y)P (y, x) si Q(x, y) = Q(y, x) = 0, et

 
π(y)Q(y, x)
π(x)P (x, y) = π(x)Q(x, y) min 1,
π(x)Q(x, y)
= min (π(x)Q(x, y), π(y)Q(y, x))
 
π(x)Q(x, y)
= π(y)Q(y, x) min ,1
π(y)Q(y, x)
= π(y)P (y, x)

si Q(x, y) et Q(y, x) sont non nuls.

Preuve de ii) : On a P (x, y) > 0 dès que Q(x, y) > 0 comme on a supposé que π(x) 6=
0. Par conséquent, toutes les trajectoires réalisables par Q le sont par P , et l’irréductibilité
de Q entraîne l’irréductibilité de P .

1
Exemple : On veut simuler sur l’espace {1, · · · n} la mesure π telle que π(k) ∝ (k+1) 2.

On se sert de la chaine de Markov auxiliaire définie comme suit : à chaque étape, lorsque
la chaine se trouve en position k, elle se déplace avec une chance sur deux sur l’un des
entiers voisins k − 1 et k + 1. Lorsque la chaine est sur un bord comme il n’y a qu’un
seul voisin disponible, avec une chance sur deux elle reste sur place, et sinon elle va sur
l’entier voisin. La matrice Q de cette chaine s’écrit alors
 
1/2 1/2
 
1/2 0 1/2 
. ...
 
Q= 1/2 . .
 

..
 

 . 0 1/2
1/2 1/2

Notons que la chaine de Markov de noyau Q est bien irréductible.

Calculons par exemple
π(2)Q(2, 1) 1 22 4
P (1, 2) = Q(1, 2) min(1, ) = min(1, 2 ) =
π(1)Q(1, 2) 2 3 18
On peut vérifier que pour tout 1 ≥ k ≤ n − 1 :
1 (k + 1)2
P (k, k + 1) =
2 (k + 2)2
pour tout 2 ≤ k ≤ n,
1
P (k, k − 1) =
2
1. ALGORITHME DE METROPOLIS-HASTINGS 51

Dans cet exemple, on remarque le fait que la matrice Q soit symétrique, simplifie
grandement les calculs.

→ Exercice 1,2,3 du TD 6 /TP 4

1.2) Espace d’état continu

Si la mesure cible π est une mesure sur R de densité f , on peut encore définir la
dynamique de type Métropolis Hastings. Dans ce cas, il faut utiliser une chaine de Markov
sur un espace non dénombrable. On se donne alors une famille de densité (q(x, y))x,y∈R
indexée par x ∈ R, c’est-à-dire que pour tout x ∈ R, on a une mesure de probabilité
q(x, y)dy, appelé noyau d’exploration. Cette famille de densités va nous permettre de
choisir un candidat pour l’évolution de l’algorithme.

Metropolis-Hastings (cas continu)

X0 = simulation selon une loi initiale quelconque ;
pour i entre 1 et n,
Y = simulation selon la loi q(Xi−1 , y)dy
f (y)q(Y,Xi−1 )
si unif orm() < f (X i−1 )q(Xi−1 ,Y )

Xi = Y
sinon faire
Xi = Xi−1
sortir (X0 , . . . , Xn );

Cas particulier : Metropolis Hastings par marche aléatoire Un cas particulier

est celui où le noyau d’exploration q(x, y) s’écrit sous la forme d’une unique loi, de densité
g : R → R, translatée par x :
q(x, y) = g(y − x).

C’est le cas dit de la marche aléatoire : si q(x, y)dy indique la loi de l’endroit où on
saute en partant de x, g désigne la loi de la taille du saut y − x.
Plus précisément, on peut remarquer que pour x fixé, la loi d’une variable Y distribuée
selon la mesure q(x, y)dy = g(y − x)dy est celle de la variable aléatoire x + Z où Z est
distribuée selon la densité g.
Plutôt que de changer de noyau d’exploration quand l’état x évolue, on peut donc toujours
utiliser la même densité g, qui nous montre la taille du prochain saut. De plus, on remarque
que dans ce cas
q(y, x) g(x − y) g(−Z)
= =
q(x, y) g(y − x) g(Z)
52 CHAPITRE V. MCMC

Metropolis-Hastings (cas de la marche aléatoire)

X0 = simulation selon une loi initiale quelconque ;
pour i entre 1 et n,
Z = simulation selon la densite g
Y = Xi−1 + Z
f (Y ) g(−Z)
si unif orm() < f (X i−1 ) g(Z)

Xi = Y
sinon faire
Xi = Xi−1
sortir (X0 , . . . , Xn );

Cas symétrique Le cas symétrique, où Q(x, y) = Q(y, x), est intéressant car certaines
simplifications apparaissent. En effet, la probabilité d’acceptation/rejet devient tout sim-
plement π(y)/π(x).
Pour le cas continu, un noyau symétrique signifie que q(x, y) = q(y, x), et donne la
probabilité d’acceptation/rejet f (y)/f (x).

Exemple : Si on veut simuler une loi de densité f (x) proportionnelle à exp(−|x|4 )

en utilisant un noyau d’exploration gaussien (on saute de x en suivant une gaussienne
standard centrée en x), on a alors un noyau d’exploration
q(x, y) = g(y − x)
où g est la densité gaussienne, qui est symétrique. On est donc dans le cas d’une marche
aléatoire symétrique et on peut écrire (en remarquant que f (y)/f (x) = exp(|x|4 − |y|4 ))
Metropolis-Hastings (cas de la marche aléatoire symétrique)
X0 = simulation selon une loi initiale ;
pour i entre 1 et n,
Z = simulation selon une gaussienne centree reduite
Y = Xi−1 + Z
si unif orm() < exp(|Xi−1 |4 − |Y |4 )
Xi = Y
sinon faire
Xi = Xi−1
sortir (X0 , . . . , Xn );

1.3) Conclusions et questions non abordées ici

L’algorithme de Metroplis Hastings est utilisé pour approcher une loi cible π qui est
connue à constante près. On a justifier sa convergence dans le cas d’un espace d’état fini.
2. ALGORITHME DU RECUIT SIMULÉ POUR LA MINIMISATION D’UNE FONCTION53

Toutefois, nous n’avons pas quantifié la vitesse de convergence de l’algorithme, qui permet
de savoir au bout de combien d’itération de la chaine de Markov, la variable simulée est
bien distribuée selon la loi invariante. Dans le cas des algorithmes de Métropolis Hastings,
cette vitesse de convergence va dépendre de la chaine de Markov auxiliaire Q choisie dans
l’algorithme.

2 Algorithme du recuit simulé pour la minimisation

d’une fonction
L’objectif de cette partie est de voir un algorithme permettant de trouver le minimum
d’une fonction définie sur un espace d’état fini (mais grand). Cet algorithme est basé
sur l’algorithme de Métropolis Hastings et la construction d’une famille de mesures de
probabilités adaptée : les mesures de Gibbs.

2.1) Mesure de Gibbs

Soit E un espace fini et V : E → R une fonction appelée traditionnellement énergie
(ou Hamiltonien) dont on cherche les points de l’espace E où le minimum est réalisé.

Définition 2.1. La mesure de Gibbs µT associée à la température T > 0 est alors la

mesure de probabilité  
1 1
µT (x) = exp − V (x)
ZT T
où ZT = y∈E exp − T1 V (y) est une constante appelée fonction de partition.
P

En général, l’espace d’états E est très grand, et la constante ZT est inconnue.

Soit M l’ensemble des points où V atteint son minimum :

M = {x ∈ E : V (x) = min V (y)}.

y

1
Proposition 2.2. i) Pour x ∈ E, lim µT (x) = .
T →∞ Card(E)
ii) Pour x ∈ E, 
1
si x ∈ M


lim µT (x) = Card(M)
T →0 
 0 sinon.

Lorsque la température T est élevée, la mesure de Gibbs a tendance à s’uniformiser

sur l’espace des états. Lorsque T diminue, µT se concentre sur les points où V atteint son
minimum.
54 CHAPITRE V. MCMC
 
1
Démonstration. 1. Pour x ∈ E, lim exp − V (x) = 1, donc
T →∞ T
  X
X 1
lim ZT = lim exp − V (x) = 1 = Card(E)
T →∞ T →∞
x∈E
T x∈E
 
1 1 1
puis lim µT (x) = lim exp − V (x) = .
T →∞ T →∞ ZT T Card(E)
2. Soit V ∗ = miny V (y). On peut écrire
∗
e−V /T
 
1 ∗
µT (x) = exp − (V (x) − V ) ,
ZT T
ce qui implique (en sommant sur les x) que
∗
e−V /T X
 
1 ∗
1= exp − (V (x) − V ) ,
ZT x∈E T

c’est-à-dire que  
ZT X 1 ∗
= exp − (V (x) − V ) .
e−V ∗ /T x∈E
T
Du coup, on peut faire lemême raisonnement que précédemment.
 Pour x ∈ E,
1 1
lim exp − (V (x) − V ∗ ) = 1 si x ∈ M et lim exp − (V (x) − V ∗ ) = 0
T →0 T T →0 T
sinon. D’où
  X
ZT X 1 ∗
lim ∗ = lim exp − (V (x) − V ) = 1 = Card(M),
T →0 e−V /T T →0 T
x∈E x∈M

ce qui permet d’écrire la limite de

∗
e−V /T
 
1 ∗
µT (x) = exp − (V (x) − V ) ,
ZT T
1
c’est-à-dire lim µT (x) = si x ∈ M et lim µT (x) = 0 sinon.
T →0 Card(M) T →0

2.2) Metropolis-Hastings pour les mesures de Gibbs

Lorsque E est grand, il peut être difficile de simuler une variable dans E selon la loi
de µT , notamment quand on ne peut pas calculer ZT . Dans ce cas-là, l’algorithme de
Metropolis-Hastings est d’un grand secours.
On considèrera dans la suite que la matrice Q est symétrique et irréductible. Dans ce
cas remarquons que  
µT (y) 1
= exp (V (x) − V (y)
µT (x) T
2. ALGORITHME DU RECUIT SIMULÉ POUR LA MINIMISATION D’UNE FONCTION55

Aussi, si V (x) > V (y) alors µµTT (x)

(y)
> 1 donc l’algorithme acceptera toujours les sauts qui
diminuent la valeur de la fonction V .
A l’inverse, si V (x) < V (y), l’algorithme pourra quand même passer de x à y. Cette
propriété évite en général à l’algorithme de rester coincer dans des minimums locaux et
permet de bien explorer tout l’espace d’état.

Metropolis-Hastings pour une mesure de Gibbs

X0 = simulation selon une loi initiale quelconque ;
pour i entre 1 et n,
y = simulation selon la loi Q(X  i−1 , y)
V (Xi−1 )−V (y)
si unif orm() < exp T
Xi = y
sinon faire
Xi = Xi−1
sortir Xn

Lorsque n est grand, la variable finale Xn a alors une distribution proche de µT .

2.3) Algorithme du recuit simulé

Nous allons maintenant voir comment modifier l’algorithme de Metropolis Hastings

précédent pour qu’il converge vers les points de E où V est minimale (on parle d’argmin de
la fonction V ). On sait que lorsque T tend vers 0, la mesure µT se concentre sur les points
où V est minimal.e L’idée du recuit simulé est alors de simuler µTn , à l’aide de l’algorithme
de Metropolis-Hastings, et de faire baisser Tn au fur et à mesure de l’algorithme.
On fixe donc une dynamique Markovienne auxiliaire Q facile à simuler, et de préférence
symétrique, puis une suite (Tn )n≥1 de températures qui tend vers 0.
Pour une matrice Q symétrique, le probabilité d’acceptation/rejet s’écrit alors

    
µT (y)Q(y, x) V (x) − V (y)
min 1, n = min 1, exp .
µTn (x)Q(x, y) Tn

On effectue donc l’algorithme suivant dans le cas symétrique.

56 CHAPITRE V. MCMC

Recuit simulé
X0 = simulation selon une loi initiale quelconque ;
pour i entre 1 et n,
y = simulation selon la loi Q(X  i−1 , y)
V (Xi−1 )−V (y)
si unif orm() < exp Ti
Xi = y
sinon faire
Xi = Xi−1
sortir Xn

Comme précédemment, dès qu’un saut se fait en direction qui diminue la fonction
V , il va être automatiquement accepté. En revanche, le recuit simulé laisse une chance à
l’algorithme de passer d’un état x à un état y même si V (y) > V (x). Toutefois, comme TN
tend vers 0, la probabilité d’effectuer un tel saut décroit avec n. On peut ainsi se sortir de
minima locaux. Les algorithmes déterministes comme par exemple la descente de gradient
risquent au contraire de se coincer dans un minimum local en refusant de ré-augmenter
la valeur de la fonction V .
Le théorème suivant, qu’on ne va pas démontrer, justifie la convergence de cet algo-
rithme.

Théorème 2.3. Pour toute fonction V , il existe une constante C qui dépend de V telle
que l’algorithme de recuit simulé avec la suite de température Tn = logCn+1 se concentre
sur l’ensemble M des points où V atteint son minimum :

lim P[Xn ∈ M] = 1.
n→∞

Calibration de l’agorithme Il faut bien se rendre compte que le choix de Tn peut

s’avérer délicat. En particulier, Tn = logCn+1 est souvent trop lent, et un mauvais choix de
la constante C peut empêcher le fonctionnement de l’agorithme (cf exo 4 en TD).
On préfère souvent prendre une suite géométrique Tn+1 = λTn avec 0.5 < λ < 1. Mais il
faut faire attention tout de même car
— si on fait baisser Tn trop vite, on accepte de moins en moins les sauts qui font
augmenter V , et on peut se retrouver coincé dans un minimum local.
— si on fait baisser Tn trop lentement,l’algorithme met trop longtemps à converger.
Chapitre VI

Approximation stochastique par

Robbins-Monro (1,5 semaines)

Pré-requis : Martingales, espérance conditionnelle

1 Cadre général
Cet algorithme a été introduit en 1950 par Robbins et Monro dans le but de rechercher
les zéros d’une fonction h. Il reste encore très utilisé. C’est un algorithme récursif ( donc
économe en mémoire), basé sur des simulations aléatoires.
Son but est d’identifier la solution θ∗ de l’équation h(θ) = 0. Il s’applique dans tous
les cas où la fonction h peut s’écrire comme une espérance

h(θ) = E(H(θ, Y ))

pour Y une variable aléatoire facile à simuler.

L’algorithme s’écrit alors

θn+1 = θn − γn+1 H(θn , Yn+1 )

avec
— (γn+1 )n≥0 une suite de pas positifs tels que
X X
γn = +∞ γn2 < ∞

En particulier la suite γn → 0. Un choix classique pour cette suite est de prendre

γn = γ1 n−α avec α ∈ (1/2, 1].
— (Yn+1 )n≥0 une suite de variables aléatoires i.i.d. de même loi que Y .

57
58 CHAPITRE VI. ROBBINS-MONRO

1.1) Exemples d’utilisation

Dans la suite, nous allons voir une série d’exemples d’utilisation, avant de démontrer
un critère de convergence de l’algorithme.

Exemple 1 : Calcul de l’espérance d’une variable aléatoire Considérons Y une

variable aléatoire, f une fonction test et la fonction

h(θ) = θ − E(f (Y )))

Résoudre l’équation h(θ) = 0 revient alors à calculer l’espérance de f (Y ). Dans ce cas, on

peut facilement écrire h comme une espérance :

h(θ) = E(θ − f (Y ))

aussi, la fonction H(θ, Y ) = θ − f (Y ).

Dans ce cas, l’algorithme de Robbins Monro s’écrit

θn+1 = θn − γn+1 (θn − f (Yn+1 )).

Nous allons voir qu’il s’agit en fait d’une généralisation de la méthode de Monte Carlo.
Prenons γn+1 = 1/(n + 1), on obtient alors

1 f (Yn+1 )
θn+1 = θn (1 − )+
n+1 n+1
soit
(n + 1)θn+1 = nθn + f (Yn+1 )

Si on pose vn = nθn , alors on a la relation de récurrence

vn+1 = vn + f (Yn+1 )

soit en fait
n
X
vn = f (Yk+1 )
k=1

ou encore
n
1X
θn = f (Yk+1 )
n k=1

Dans ce cas, θn est la moyenne empirique (ou estimateur de Monte-Carlo).

1. CADRE GÉNÉRAL 59

Exemple 2 : Estimation d’un quantile On rappelle que si Y a pour fonction de

répartition F , on appelle le quantile d’ordre α la valeur θα telle que F (θα ) = α.
En particulier, le quantile est solution de

0 = F (θ) − α
= E(1Y ≤θ − α)

Pour écrire l’agorithme de Robbins Monro, on pose H(θ, Y ) = 1Y ≤θ − α et on obtient

θn+1 = θn − γn+1 (1Yn+1 ≤θn − α)

L’avantage de cet algorithme par rapport aux autres méthodes d’estimation de quan-
tile, notamment basées sur des méthodes de tri est son coup de calcul très faible.

Descente de gradient stochastique On peut utiliser l’algorithme de Robbins Monro

pour trouver les minimas (ou maximas) d’une certaine fonction.
C’est le cas typiquement lorsque la fonction h correspond à un gradient (ou à l’opposé
d’un gradient) :
h(θ) = −∇V (θ) = E[H(Y, θ)]
(en dimension 1, ∇V est juste la dérivée de V ).
Un minimiseur de V est alors aussi une solution de h(θ) = 0, et on peut utiliser l’algo-
rithme de Robbins-Monro pour le trouver.
Cela conduit à des modifications stochastiques de l’algorithme de descente de gradient
déterministe.
- Exemple (gradient bruité) : c’est le cas où on connaît ∇V (θ) à un aléa près, par
exemple si H(Y, θ) = −∇V (θ) − Y avec E[Y ] = 0. On est bien dans le cas où

E[H(Y, θ)] = h(θ) = −∇V (θ),

et la récurrence de Robbins-Monroe s’écrit alors

θn+1 = θn − γn+1 (∇V (θn ) + Yn+1 ),

où (Yn )n≥1 est une suite de copies de Y indépendantes et (γn )n≥0 une suite de pas qui
tend vers 0.
→ En TD on appliquera ceci pour la regression linéaire.

Exemple d’utilisation dans un cadre non stationnaire : Dans les exemples pré-
cédents, les variables aléatoires (Yn )n≥1 étaient identiquement distribuées.
Nous allons passer à une situation où nous ne sommes plus seulement observateur
de la quantité aléatoire (Yn )n≥1 , mais aussi acteur (imaginez une température intérieure
60 CHAPITRE VI. ROBBINS-MONRO

qui dépend de la température de vos radiateurs, ou un taux d’une certaine quantité

physiologique contrôlée par l’injection d’un médicament). Le paramètre qui nous permet
d’agir sur les (Yn )n≥1 est justement le paramètre θ, qu’on cherche à ajuster.
On se donne donc une famille de loi (pθ )θ et on va considérer une variable aléatoire Y
de loi pθ (il faut donc mettre l’indice θ à chaque fois qu’on calcule une probabilité ou une
espérance sur Y ).
On cherche comme précédemment la solution θ∗ de h(θ) = 0, où h(θ) = Eθ [H(Y, θ)].
On va donc considérer la suite d’estimateur de θ∗ définie par récurrence par

θn+1 = θn + γn+1 H(Yn+1 , θn ),

où (Yn )n≥1 est une suite de variables aléatoires de lois conditionnelles fixées : pθn
est la loi conditionnelle Yn+1 sachant θn ; et (γn )n≥0 est une suite de pas qui tend vers 0.
On peut l’écrire comme ceci.

Algorithme de Robbins-Monro (non-stationnaire

θ0 = arbitraire ;
pour n entre 0 et N − 1,
Yn+1 = simulation selon la loi pθn
θn+1 = θn + γn+1 H(Yn+1 , θn )
sortir θN
Sous des hypothèses que l’on ne précisera pas ici,

lim θn = θ∗ ps
n→∞

et
θn − θ∗
lim √ = gaussienne en loi.
n→∞ γn

Exemple (dosage) : Y est une grandeur physiologique observable chez un patient

auquel on injecte une quantité de produit θ. On cherche la quantité idéale θ∗ de produit
à administrer au patient pour obtenir l’effet moyen escompté α sur Y . Autrement dit, on
cherche la solution θ∗ de
Eθ [Y ] = α.
En posant H(Y, θ) = α − Y , on obtient

h(θ) = E[H(Y, θ)] = E[α − Y ] = α − E[Y ],

qui s’annule en θ∗ . La récurrence de Robbins-Monro s’écrit alors

θn+1 = θn + γn+1 (α − Yn+1 ),

2. CONVERGENCE DE L’ALGORITHME 61

où sachant θn , Yn+1 est tirée sous la loi pθn .

Si Eθ [Y ] est croissant en θ (action positive de θ sur Y ), la situation est assez intuitive.
Si Y est trop bas, on augmente θ et si Y est trop haut, on diminue θ.
On verra un second exemple en TP pour la régression linéaire sur un modèle de données
auto-régressives.

2 Convergence de l’algorithme
Pour rappel, cherche à trouver la solution θ∗ de l’équation h(θ) = 0, avec θ un réel (ou
à valeurs dans Rd ) mais on n’a pas accès à h directement. On sait que

h(θ) = E[H(Y, θ)],

où Y est une variable aléatoire et θ 7→ H(Y, θ) une version aléatoire de h. On a accès à

des réalisations indépendantes de variables aléatoires distribuées comme Y .
L’algorithme de Robbins Monro calcule la suite définie par récurrence par

θn+1 = θn + γn+1 H(Yn+1 , θn ),

où (Yn )n≥1 est une suite de copies de Y indépendantes et (γn )n≥0 une suite de pas qui
tend vers 0.

P P 2
Si n γn = +∞ et n γn < +∞, et sous de bonnes hypothèses sur H, on a les
convergences suivantes :
lim θn = θ∗ ps
n→∞
et
θn − θ∗
lim √ = gaussienne en loi.
n→∞ γn
Dans cette partie du cours nous allons donner un cadre pour la première convergence
et nous illustrerons en TP la seconde convergence.
Pour démonter la convergence presque sure de l’algorithme de Robbins-Monro, l’idée
principale est d’utiliser la structure récursive de cet algorithme et de mettre en évidence
l’existence d’une sur-martingale positive associée.
Pour cela, on a besoin de définir la filtration associée à l’algorithme. Pour tout n ≥ 0,
on définit Fn la tribu engendrée par les variables aléatoires (θ0 , · · · , θn , Y0 , · · · Yn ).
On remarque alors que

E(θn+1 |Fn ) = E(θn − γn+1 H(θn , Yn+1 )|Fn )

= θn − γn+1 E(H(θn , Yn+1 )|Fn )
= θn − γn+1 h(θn )
62 CHAPITRE VI. ROBBINS-MONRO

car θn est Fn mesurable est Yn+1 est indépendante de la tribu Fn . On peut donc écrire

θn+1 = θn − γn+1 h(θn ) + γn+1 ∆Mn+1

avec
∆Mn+1 = h(θn ) − H(θn , Yn+1 ))

On appelle la quantité ∆Mn+1 un incrément de martingale car E(∆Mn+1 |Fn ) = 0. En

particulier Mn = nk=1 ∆Mk sera une martingale.
P

L’algorithme fait évoluer la valeur de θn en suivant un champ aléatoire θ 7→ H(Y, θ)

et on appelle naturellement champ moyen la fonction h.
Lorsque les pas sont petits, l’algorithme a tendance à suivre les trajectoires moyennes,
c’est à dire les solutions θ̄ : R+ → R de l’équation différentielle θ̄0 = h(θ̄) .
La convergence de l’algorithme est en fait une conséquence de la convergence des
trajectoires moyennes vers les solutions de h(θ) = 0 pour tout point de départ θ̄(0).

2.1) Théorème de Robbins-Siegmund

Nous allons commencer par introduire un résultat général très puissant pour l’étude
des algorithmes stochastiques.

Théorème 2.1 (Théorème de Robbins-Siegmund). Soient (Un ), (Vn ), (αn ) et (βn ) quatre
suites de variables aléatoires positives et Fn -adaptées telles que
(1) (αn ), (βn ) et (Un ) sont prévisibles (i.e. αn+1 est Fn − mesurable).
(2) supn∈N nk=1 (1 + αk ) < ∞ presque surement
Q
P∞
(3) n=0 E(βn ) < ∞
(4) Pour tout n ∈ N,

E(Vn+1 |Fn ) ≤ Vn (1 + αn+1 ) + βn+1 − Un+1 (VI.1)

Alors
a) Vn −→ V∞ presque surement et V∞ ∈ L1
b) supn∈N E(Vn ) < ∞
P P
c) n E(Un ) < ∞ et n Un < ∞ presque sûrement.

Le reste de la section est consacrée à la preuve de ce résultat. C’est une preuve tech-
nique mais instructive, car elle détaille la construction d’une sur-martingale positive et
exploite le fait qu’une sur-martingale positive converge presque surement vers une variable
aléatoire intégrable.

Démonstration. Construction de la sur-martingale :

2. CONVERGENCE DE L’ALGORITHME 63
Pn
On considère tout d’abord la quantité (Vn + k=1 Uk ). En utilisant (4), on obtient

n+1
! n+1
X X
E Vn+1 + Uk |Fn = E (Vn+1 |Fn ) + Uk
k=1
|k=1{z }
prévisible
n+1
X
≤ Vn (1 + αn+1 ) + βn+1 − Un+1 + Uk
k=1
n
X
≤ Vn (1 + αn+1 ) + βn+1 + Uk
k=1
n
X
≤ (Vn + Uk )(1 + αn+1 ) + βn+1
k=1

Définissons maintenant Pn
Vn + Uk
k=1
Sn = Qn (VI.2)
k=1 (1 + αk )

alors
βn+1
E(Sn+1 |Fn ) ≤ Sn + Qn+1 (VI.3)
k=1 (1 + αk )

Considérons la suite
n
X βk
Bn = Qk
k=1 j=1 (1 + αj )

La variable Bn est une somme de terme positifs, aussi la suite Bn est croissante et elle
converge presque surement vers B∞ (une quantité qui peut être infinie). On veut montrer
que E(B∞ ) est finie. Pour ce faire, remarquons que

E(B∞ ) = E(lim Bn ) = lim E(Bn )

n n

par théorème de convergence monotone. De plus

n
! n ∞
X βk X X
E(Bn ) = E Qk ≤ E(βk ) ≤ E(βk ) < ∞
k=1 j=1 (1 + αj ) k=1 k=1

par l’hypothèse (3). Donc E(B∞ ) < ∞ et par conséquent la variable B∞ est finie presque
surement.
Finalement on définit
Sen = Sn + E(B∞ |Fn ) − Bn .
64 CHAPITRE VI. ROBBINS-MONRO

Remarquons d’abord que Sen ≥ 0 car Sn ≥ 0 et E(B∞ |Fn ) − Bn = E(B∞ − Bn |Fn ) ≥ 0

puisque Bn est croissante vers B∞ . De plus, avec les inégalités précédentes :

E(Sen+1 |Fn ) = E(Sn+1 + E(B∞ |Fn+1 ) − Bn+1 |Fn )

βn+1
≤ Sn + Qn+1 + E(B∞ |Fn ) − Bn+1
k=1 (1 + α k )
≤ Sn + E(B∞ |Fn ) − Bn = Sen

On en déduit donc que Sen est une sur-martingale positive et donc qu’elle converge vers
une variable aléatoire Se∞ qui est intégrable.
Conclusions : Ecrivons maintenant

Sn = Sen − (E(B∞ |Fn ) − Bn )

d’une part Sen → Se∞ presque surement et par croissance (E(B∞ |Fn ) − Bn ) → 0 presque
surement.
Ainsi
Sn → Se∞ .

On veut maintenant revenir à Vn et à nk=1 Uk pour démontrer le théorème.

P

En prenant l’espérance dans (VI.2) on peut écrire

Xn n
Y
E(Vn ) + E( Uk ) = E(Sn (1 + αj )).
k=1 j=1

Aussi chacun des termes du membre de gauche est inférieur au membre de droite. Or
n
Y n
Y
sup E(Sn (1 + αj )) ≤ sup (1 + αj )) sup E(Sn )
n n n
j=1 j=1
| {z }
<∞p.s. (2)

Or
E(Sn ) ≤ E(Sen ) ≤ E(Se1 ) < ∞

Donc simultanément
Xn
sup E(Vn ) < ∞ sup E( Uk ) < ∞
n n
k=1

Comme les Un sont positives on en déduit que

X∞ ∞
X
E( Uk ) < ∞ et Uk < ∞ .p.s.
k=1 k=1
2. CONVERGENCE DE L’ALGORITHME 65

Ce qui démontre les points b) et c) du théorème.

Finalement n n
Y X
Vn = Sn (1 + αj ) − Uk
j=1 k=1

comme tous les termes du membre de droite convergent presque surement vers une variable
L1 on en déduit que Vn −→ V∞ et E(|V∞ |) < ∞.

2.2) Application à Robbins Monro

Dans cette section, on va voir comment ce résultat nous permet de démontrer la
convergence presque sûre de l’algorithme de Robbins-Monro. On rappelle que cet algo-
rithme s’écrit
θn+1 = θn − γn+1 h(θn ) + γn+1 ∆Mn+1
avec l’incrément de martingale ∆Mn1+ = h(θn ) − H(θn , Yn+1 ) qui vérifie E(∆Mn+1 |Fn ) =
0.
La preuve fait intervenir une fonction de Lyapunov pour l’algorithme

Définition 2.2. Dans ce cours, on dira qu’une fonction V : E 7→ R+ est une fonction de
Lyapunov si
i) pour tout (x, y) ∈ E 2 , ||∇V (x) − ∇V (y)||≤ L||x − y|| pour L > 0 une constante.
(gradient -Lipschitz)
ii) minx∈E V (x) = m > 0
iii) lim|x|→∞ V (x) = +∞ (coercive)
iv) ||∇V (x)2 ||≤ C(1 + V (x)) (sous quadratique)

Un exemple typique de telle fonction est la fonction V (x) = x2 + 1.

Nous allons maintenant donner les conditions de convergence presque sure de l’algo-
rithme.

Théorème 2.3. Supposons qu’il existe une fonction de Lyapunov V vérifiant les hypo-
thèses de la Définition 2.2 et se comportant bien vis à vis de l’algorithme, c’est à dire telle
que :
1) h∇V (x), h(x)i ≥ 0 pour tout x ∈ E (positive le long des trajectoires)
2) |h(x)|2 ≤ C(1 + V (x)) (asymptotique contrôlée par V )
3) E(|∆Mn+1 |2 |Fn ) ≤ C(1 + V (θn )) (contrôle du bruit)
Alors
a) supn E(V (θn )) < ∞
66 CHAPITRE VI. ROBBINS-MONRO
P
b) γn h∇V (θn ), h(θn )i < ∞
c) V (θn ) → V∞ presque surement et V∞ ∈ L1
d) θn+1 − θn → 0 presque surement et dans L2

Proposition 2.4. Si de plus

{x, h∇V (x), h(x)i = 0} = {x∗ }

alors
e) x∗ est l’unique minimum de V
7 x∗ presque surement.
f ) θn →

Preuve du Théoreme 2.3. L’idée principale de la preuve est d’effectuer un développement

de Taylor à l’ordre 2, de V (θn+1 ) au voisinage de θn est d’utiliser les hypothèses pour
obtenir une inégalité du type (VI.1).
En utilisant un développement de Taylor on a

V (θn+1 ) = V (θn − γn+1 (h(θn ) − ∆Mn+1 ))

1
= V (θn ) − γn+1 h∇V (θn ), (h(θn ) − ∆Mn+1 )i + D2 V (ξn+1 ).(γn+1 (h(θn ) − ∆Mn+1 )))2
2

Le contrôle du terme d’ordre 2 provient du fait que V est gradient Lipchitz, donc sa
hessienne est bornée :
1 2
D V (ξn+1 ).(γn+1 (h(θn ) − ∆Mn+1 )))2 ≤ Cγn+1
2
(h(θn ) − ∆Mn+1 ))2
2
2
≤ Cγn+1 (h(θn )2 + ∆Mn+1 )2 )

Donc en utilisant 2) et 3),

1
E D2 V (ξn+1 ).(γn+1 (h(θn ) − ∆Mn+1 )))2 |Fn ) ≤ Cγn+1
2
(1 + V (θn ))
2
Pour les termes d’ordre 1, remarquons que comme ∆Mn+1 est un incrément de martingale

E(h∇V (θn ), (h(θn )−∆Mn+1 )i|Fn ) = h∇V (θn ), (h(θn )−E(∆Mn+1 )|Fn )i = h∇V (θn ), h(θn )i

D’où finalement
2
E(V (θn+1 )|Fn ) ≤ V (θn ) − γn+1 h∇V (θn ), h(θn )i + Cγn+1 (1 + V (θn ))
2 2
≤ V (θn )(1 + Cγn+1 ) − γn+1 h∇V (θn ), h(θn )i + Cγn+1 (VI.4)

On peut donc appliquer le théorème de Robbins Siegmund avec Vn = V (θn ), an+1 =

2
cγn+1 = βn+1 et Un+1 = γn+1 h∇V (θn ), h(θn )i. Vérifions les hypothèses du théorème :
toutes ces suites sont positives, adaptées et prévisibles quand il le faut 1). Par hypothèse
P 2
sur γn+1 < ∞, on en déduit que 2) et 3). On obtient alors les conclusions suivantes :
2. CONVERGENCE DE L’ALGORITHME 67

— V (Xn ) 7→ V∞ , et V∞ ∈ L1 ce qui donne c)

P
— Un < ∞ ce qui donne b)
— supn E(V (θn )) < ∞ ce qui donne a).
Il reste à montrer d). On sait par hypothèse que

E((θn+1 − θn )2 |Fn ) ≤ Cγn+1

2
(1 + V (θn ))

donc X
E((θn+1 − θn )2 ) < ∞.
On en déduit que E((θn+1 − θn )2 ) → 0 mais aussi que (θn+1 − θn )2 < ∞ presque
P

surement, donc θn+1 − θn → 0 presque surement et dans L2 .

Preuve de la proposition 2.4. Pour le premier point remarquons déjà que tous les minimas
de V sont inclus dans l’ensemble {x, h∇V (x), h(x)i = 0}. Donc si cet ensemble est un
singleton, V a un unique minimum.
On sait maintenant que V atteint son minimum en x∗ . Montrons que lim V (θn ) = V∞ =
V (x∗ ) et lim ∇V (θn ) = 0 on aura alors lim θn ∈ {x, h∇V (x), h(x)i = 0} ce qui conclura
la preuve.
On se place sur l’événement Ω de probabilité 1, sur lequel V (θn (ω)) → V∞ (ω) et
X
γn+1 h∇V (θn (ω)), h(θn (ω))i < ∞.

En particulier, comme V (θn (ω)) converge, la suite θn (ω) est bornée, donc elle converge
à extraction près vers θ∞ (ω). On déduit de la convergence de la série, que nécessairement

limh∇V (θn (ω)), h(θn (ω))i = 0,

et donc par continuité que

h∇V (θ∞ (ω)), h(θ∞ (ω))i.
Par hypothèse, θ∞ (ω) = x∗ . En particulier, θn (ω) a une unique valeur d’adhérence et
converge donc vers cette valeur. Finalement, V∞ = V (x∗ ).

On verra un exemple d’utilisation de ce résultat en TD (Exercice 1 sur les quantiles).

Nous avons justifié de la convergence de l’algorithme. Nous verrons en TP deux types

de résultats pour préciser cette convergence :
— un TCL sur la convergence en loi des fluctuations
θn − θ∗
lim √ = loi gaussienne
n→∞ γn
— des estimées sur le comportement à n fixé, de l’erreur moyenne quadratique

E((θn − θ∗ )2 ).