République Algérienne Démocratique et Populaire

Ministère de l'Enseignement Supérieur et de la Recherche

Scientique
Université Abou Bakr Belkaid- Tlemcen

Faculté des Sciences

Département des Mathématiques

MEMOIRE DE MASTER
EN MATHEMATIQUES

Option : Probabilités-Statistiques

Méthodes de Monte-Carlo

Présenté par : MEDOUAKH Fatima Zahra

Mémoire soutenu devant le jury composé de :

M. T. Mourid Professeur Universsité de Tlemcen Président

Mme. M. Dali youcef M.C.A Universsité de Tlemcen Encadreure
M. A. Allam M.C.A Universsité de Tlemcen Examinateur
M. A. Labbas M.C.A Universsité de Tlemcen Examinateur

Année universitaire 2016-2017

 Remerciements

Tout d'abord, je remercie ALLAH qui m'a donné la volonté, la patience,

et surtout la santé durant toutes mes années d'étude.

Mes vifs remerciements sont adressés à mon encadreure Malika DALI

YOUCEF pour sa patience, ses conseils, ses orientations, et sa disponibilité
durant la préparation de ce mémoire de n d'études de Master.

Je voudrais remercier tous les professeurs qui ont contribué à ma

formation et à la réussite de mes études universitaires ainsi que les membres

de mon Jury de soutenance :

monsieur Tahar MOURID qui m'a fait le grand honneur d'accepter, d'être
le président de mon jury, monsieur Ahmed LABBAS et monsieur Abdelaziz
ALLAM qui ont pris la peine de me lire et d'évaluer mon travail.

Enn ma pensée et ma gratitude vont à l'ensemble des responsables et

de l'administration du Département de Mathématiques, sans oublier tout
le personnel employé au Département .

Je tiens aussi à remercier tout particulièrement mes parents d'abord et

ma famille ensuite pour leur encouragements continus.

i
Table des matières

1 Introduction iii
2 Simulation de variables aléatoires. 1
2.1 Simulation de la loi uniforme . . . . . . . . . . . . . . . . . 1
2.2 Méthode d'inversion . . . . . . . . . . . . . . . . . . . . . . . 1
2.3 Méthode de rejet . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Simulation des lois normales . . . . . . . . . . . . . . . . . . 10

3 Calcul d'une intégrale par la méthode de Monte-Carlo. 13

4 Méthode de Monte-Carlo par Chaînes de Markov. 20
4.1 La propriété de Markov . . . . . . . . . . . . . . . . . . . . . 20
4.2 Calcul des lois marginales . . . . . . . . . . . . . . . . . . . . 21
4.3 Récurrence et transience . . . . . . . . . . . . . . . . . . . . . 24
4.4 Mesure invariante . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.5 Période d'une chaîne de Markov . . . . . . . . . . . . . . . . . 30
4.6 Méthode de Monte-Carlo . . . . . . . . . . . . . . . . . . . . 31

5 Conclusion 33
6 Annexe 34

ii
1 Introduction

La méthode de Monte-Carlo a vu son essor à partir de la n de la seconde

guerre mondiale, essentiellement dans le cadre du projet americain "Manhat-
tan" concernant le développement de l'arme nucléaire. Cette époque corres-
pond également à la construction des premiers "ordinateurs". Ce projet étant
été classé "secret défense", il est dicile de savoir exactement qui parmi ses
pionniers : Von Neumann, Ulam, Metropolis a proposé le nom de "Monte-
Carlo". Quoi qu'il en soit, ce terme fait référence aux jeux de hasard : la
capitale de la principauté de Monaco.
Le terme méthode de Monte-Carlo désigne une famille de méthodes algo-
rithmiques visant à calculer une valeur numérique approchée en utilisant
des procédés aléatoires, c'est-à-dire des techniques probabilistes. D'un point
de vue purement mathématique, une méthode de Monte-Carlo peut servir
au calcul d'intégrales (simples ou multiples) et à la résolution d'équations
aux dérivées partielles, de systèmes linéaires et de problèmes d'optimisation.
D'une manière plus globale, elle permet d'avancer toutes les possibilités liées
à une observation, en quantiant le risque associé (ou l'incertitude statis-
tique). Pratiquement, ces techniques sont couramment utilisées dans divers
domaines (physique, ingénierie, nance. . . ).
Dans ce mémoire, nous avons présenté la méthode de Monte-Carlo du point
de vue de son application pour le calcul d'intégrale. En particulier, pour ap-
pliquer la méthode en question, il faut savoir simuler les variables aléatoires,
c'est la raison pour la quelle nous avons introduit dans la première section
la simulation de variables aléatoires dont la méthode d'inversion, méthode
de rejet, ainsi que la méthode de Box-Muller.
Le problème traité est donc le calcul de h(x)dπ(x), où π est une loi de
R

probabilité et h est une fonction donnée. Il existe des cas où on ne peut pas
simuler π par les méthodes précédentes de simulation, la méthode de Monte
Carlo par chaîne de Markov pour résoudre de tels problèmes "plus nis" et
c'est le thème de la troisième section.

iii
2 Simulation de variables aléatoires.

Pour appliquer une méthode de Monte-Carlo, il faut savoir simuler sui-

vant une loi donnée, dans cette section nous allons voir quelques méthodes
de simulation de variables aléatoires.

2.1 Simulation de la loi uniforme

La distribution de base qu'il faut savoir simuler numériquement est la

loi uniforme sur [0, 1]. Tous les langages de programmation disposent d'un
générateur de nombres au hasard dans [0, 1]. Par exemple, en Matlab le gé-
nérateur est rand.
Maintenant, pour simuler X une variable aléatoire qui suit la loi uniforme
sur [a,b], il sut de simuler Y = a + (b − a)U , où U suit la loi uniforme sur
[0, 1].
Il existe plusieurs méthodes de simulation : Nous rappelons dans le para-
graphe suivant quelques une de ces méthodes de simulation.

2.2 Méthode d'inversion

principe

La méthode d'inversion est la plus simple des méthodes générales de

simulation. Elle consiste à composer l'inverse de la fonction de répartition F
de la distribution à simuler avec un générateur de la loi uniforme sur [0,1].
F n'étant pas toujours inversible au sens classique, on déni son inverse de
la façon suivante :

∀u ∈ [0, 1], F −1 (u) = inf{x; F (x) ≥ u}.

De point de vue théorique, le principe de la méthode est justié dans la

proposition suivante :
Proposition 2.1. Soient F une fonction de répartition sur R et U une
variable aléatoire de loi uniforme sur [0, 1]. Alors la variable aléatoire
X = F −1 (U ),

a pour fonction de répartition F .

1
 Pour la preuve nous avons besoin du lemme suivant :
Lemme 2.1. Pour tout u ∈ [0, 1], x ∈ R
F −1 (u) ≤ x ⇐⇒ u ≤ F (x).

Démonstration du lemme. Pour l'implication directe, la croissance de F

donne
F (F −1 (u)) ≤ F (x),
puisque la fonction de répartition est continue à droite
F (F −1 (u)) ≥ u,

nous avons donc F (x) ≥ u, quant à l'implication réciproque elle est triviale.

Démonstration de la proposition.
∀x ∈ R, P (X ≤ x) = P (F −1 (U ) ≤ x) = P (U ≤ F (x)) = F (x).

Cette proposition permet d'obtenir les formules analytiques utiles pour

simuler un grand nombre de lois, nous introduisons quelques exemples illus-
tratifs :

Exemples

1. La loi exponentielle : soit X une variable aléatoire réelle qui suit une
loi exponentielle de paramètre λ. La fonction réciproque de sa fonction
de répartition vaut pour u ∈ [0, 1],
− ln(1 − u)
F −1 (u) = ,
λ
et si U est une variable aléatoire de loi uniforme sur [0, 1], X est de
même loi que
− ln(1 − U )
,
λ
ou encore
− ln U
.
λ
Le code R suivant donne 100 observations de la loi exponentielle de
paramètre 1

2
 N = 100
U =runif(N )
X = −log(1 − U )

Figure 1  100 observations de la loi exponentielle de paramètre 1.

2. La Loi de Cauchy : soit X une variable aléatoire réelle de loi de Cauchy

de paramètre a, alors la variable
1
F −1 (U ) = a tan(π(U − )),
2
suit une loi de Cauchy de paramètre a.
Nous avons proposé le programme suivant qui permet de simuler 100
observations de la loi de Cauchy de paramètre 1 :
N = 100
U = runif (100)
X = numeric(N )
f or(i in 1 : N )
{X[i] = tan(pi ∗ (U [i] − 1/2))}
plot(X)

3
 Figure 2  100 observations de la loi de Cauchy de paramètre 1.

On peut aussi simuler les variables discrètes par la méthode d'inversion.

Lois discrètes

Soit X une variable aléatoire réelle discrète de loi pi = P (X = xi )., nous

supposons que les éventualités sont rangé par ordre croissant

xi < xi+1 , i ∈ {1...n}, n ∈ N,

considérons La fonction de répartition de X est dénie par :

0 si x < x1

F (x) =
p1 + ...pi = Fi si xi ≤ x < xi+1

Nous pouvons dénir l'algorithme de simulation par la méthode d'inversion

comme suit :
i←1
choix ← Random
Tant que (choix > Fi ) faire
i←i+1
n Tant que
X ← xi .

 Exemple : simulation de la loi de Poisson

4
Si X suit une loi de Poisson de paramètre λ, il n'y a pas d'expression simple
pour la fonction de répartition et l'ensemble des valeurs possibles est inni.
Il faut donc calculer les valeur Fi au fur et à mesure.
L'algorithme relatif à une telle simulation est donc :
X←0
P ← exp −λ
F ←P
choix ← Random
Tant que (choix > F ) faire
X ←X +1
P ← P λ/X
F ←F +P
n Tant que
X ← X − 1.

Remarquons que la méthode d'inversion n'est exacte qu'à condition de

connaitre l'expression explicite de F −1 . C'est rarement la plus ecace pour
les variables à densité. Une seconde méthode de simulation est introduite
dans le paragraphe suivant.

2.3 Méthode de rejet

La méthode de rejet est utilisée pour engendrer indirectement une va-

riable aléatoire X , de densité de probabilité f , lorsqu'on ne peut pas utiliser
la simulation par inversion. Le plus simple est de l'introduire en premier pour
une loi uniforme.

Méthode de rejet pour la loi uniforme

On suppose que l'on sait simuler une loi uniforme surD, et que l'on veut
simuler une loi uniforme sur C ⊂ D. Alors l'algorithme de simulation par la
méthode de rejet est le suivant :
répéter
tirer un point X au hasard dans D
jusqu'à (X ∈ C).

Pour le choix du support D, on a intérêt à le choisir le plus proche possible

de C, mais susamment simple pour ne pas ralentir la simulation. Citons
comme exemple la simulation de la loi uniforme sur le disque unité :
C = le disque unité.

5
On choisit D = le carré [−1, 1]2 .
L'algorithme est :
répéter
X ← 2 Random − 1
Y ← 2 Random − 1
S ← X2 + Y 2
jusqu'à (S ≤ 1).

La simulation des lois uniformes revient souvent dans les applications. En

ce sens que le résultat suivant montre que la simulation d'une loi de densité
quelconque peut toujours se ramener à la simulation d'une loi uniforme.
Proposition 2.2. Soit f une densité de probabilité par rapport à la mesure
de Lebesgue sur Rd , continue par morceaux, et
D = {(x, u) ∈ Rd × R+ ; 0 < u < f (x)}.

Soit X un vecteur aléatoire à valeurs dans Rd et U une variable aléatoire

réelle. le couple (X,U) suit alors la loi uniforme sur D si et seulement si :
1. X a pour densité f.

2. La loi conditionnelle de U sachant "X = x" est la loi uniforme sur

[0, f (x)].

Démonstration. Pour l'implication direct , (X, U ) suit une loi uniforme sur
D, alors sa fonction de densité est :
1
f(X,U ) (x, u) = 1D (x, u),
vol(D)

vol(D) = 1 car f est une densité.

Pour montrer le point 1 de la proposition, on intègre f(X,U ) par rapport à u

Z Z
fX (x) = f(X,U ) (x, u)du = 1D (x, u)du
R R
Z f (x)
= du = f (x),
0

6
pour le point 2 :
f(X,U ) (x, u)
fUX=x (u) =
fX (x)
1
= 1D (x, u)
f (x)
1
= 1 (u).
f (x) [0,f (x)]

Et pour l'implication inverse, remarquons que la loi du couple (X, U ) est le

produit de la loi de X avec la loi conditionnelle de U sachant "X = x"

f(X,U ) (x, u) = fX (x)fUX=x (u)

1[0,f (x)] (u)
= f (x)
f (x)
= 1[0,f (x)] (u).

En d'autres termes, une variable aléatoire de densité f est l'abscisse d'un

point au hasard sous le graphe de f . Ce résultat établi pour des loi à densités
par rapport à la mesure de Lebesgue, s'étend naturellement à des mesures
quelconques.
Dans le paragraphe suivant, la méthode de rejet est appliquée pour la simu-
lation de loi à densité.

Méthode de rejet pour les lois à densité

Nous allons voir une proposition qui permet à partir d'une densité g facile
à simuler, simuler une loi de densité f.
Proposition 2.3. Soit µ une mesure positive sur Rd , f et g deux densités
de probabilité sur (Rd , µ) telles qu'il existe une constante c vériant :

∀x ∈ Rd , cg(x) ≥ f (x).
Soit X une variable aléatoire de densité g (par rapport à µ) et U une variable
aléatoire de loi uniforme sur[0, 1], indépendante de X . Alors la loi condition-
nelle de X sachant l'événement {cU g(X) < f (X)}, a pour densité f par
rapport à µ.

7
Démonstration. On cherche la loi conditionnelle de X sachant l'évènement
A = {cU g(X) < f (X)},

pour cela calculons d'abord la probabilité de A.

Remarquons que c ≥ 1 car
Z Z
c g(x)dµ(x) ≥ f (x)dµ(x) = 1,
Rd Rd

donc c 6= 0 et cg(x)
f (x)
≤ 1.
D'autre part, le support de f est inclus dans le support de g , alors

P (A) = P ((X, U ) ∈ {(x, u), cug(x) < f (x)})

Z
= g(x)1[0,1] (u)dudµ(x)
{(x,u),cug(x)<f (x)}
Z Z f (x)
cg(x)
= g(x)( du)dµ(x)
{x,g(x)>0} 0
Z
f (x) 1
= g(x) dµ(x) = .
{x,g(x)>0} cg(x) c

On veut montrer que la densité conditionnelle de X sachant A est f , c'est-

à-dire que pour tout borélien B de Rd :
Z
P (X ∈ B/A) = f (x)dµ(x)
B

P (X ∈ B/A) = c P ({X ∈ B} ∩ A)
Z Z f (x)
cg(x)
=c g(x)( du)dµ(x)
{x,g(x)>0}∩B 0

par le calcul précédent, nous avons

Z Z
P (X ∈ B/A) = f (x)dµ(x) = f (x)dµ(x)
{x,g(x)>0}∩B Rd ∩B
Z
= f (x)dµ(x).
B

8
 L'algorithme de simulation de la méthode de rejet est alors :
répéter
simuler X de densité g
U ← Random
jusqu'à (cU g(X) < f (X)).

On remarque que l'algorithme comporte une boucle dont la condition

porte sur des variables aléatoires. le nombre d'itérations, noté N est donc lui-
même aléatoire. On peut montrer que N suit la loi géométrique de paramètre
c , car :
1

P (N = k) = (1 − p)k−1 p,
d'après la preuve de la proposition précédente nous avons :
1
p = P (cU g(X) < f (X)) = .
c
L'espérance de N c'est-à-dire le nombre moyen d'itération à eectuer avant
d'obtenir une réalisation de la densité f vaut c.
E(N ) = c.

On a donc tout intéret à choisir c le plus petit possible. Comme exemple

de la simulation par la méthode de rejet, citons la simulation de la densité
gaussienne centrée réduite à partir de la densité de Laplace.
Tout d'abord pour la simulation de la densité de Laplace, nous pouvons
utiliser la méthode d'inversion :
La densité de Laplace est :
e−|x|
g(x) = , x ∈ R,
2
et sa fonction de répartition :
ex
si x < 0

G(x) = 2
e−x
1− 2 si x ≥ 0
son inverse est déni comme suit :
ln(2u) si 0 < u < 1

−1 2
G (u) =
− ln(2(1 − u)) si u ≥ 12

L'algorithme de simulation de g est le suivant :

U ← Random

9
si U < 1
2
X ← ln(2U )
sinon
X ← − ln(2(1 − U )).

Soit f la densité gaussienne, nous avons :

r
2e
f (x) ≤ g(x), ∀x ∈ R.
π
En eet : 1 2 1 1 1
2
e− 2 x +|x|
= e 2 e− 2 (|x|−1) ≤ e 2 ,
r
1 1 2 2e e−|x|
=⇒ √ e− 2 x ≤ ,
2π π 2
q
alors il existe c = 2e
π tel que :

f (x) ≤ c g(x),

nous pouvons alors appliquer l'algorithme de simulation par rejet.

2.4 Simulation des lois normales

Pour la simulation de lois gaussiennes, nous rappelons l'algorithmes de

Box-Muller consiste à engendrer des couples de variables indépendantes en
simulant leurs coordonnées polaires. Le principe de la base est établi dans la
proposition suivante :
Proposition 2.4. Soit (X,Y) un couple de variables aléatoires indépen-
dantes, de même loi N (0, 1), et (R, Θ) le couple de coordonnées polaires
correspondant :
X = R cos Θ; Y = R sinΘ.
Les variables aléatoires R et Θ sont indépendantes, le module R a pour den-
sité :
fR (r) = r exp(−r2 /2)1R+ (r).
et l'argument Θ suit la loi uniforme sur [0, 2π].
Démonstration. Posons :
φ : R2 −→ R+ × [0, 2π]

10
 (x, y) −→ (r, θ)
r et θ les coordonnées polaires associées à (x,y), alors φ est bijective.

∀(r, θ) ∈ R+ × [0, 2π], φ−1 (r, θ) = (r cos θ, r sin θ),

Jac(φ−1 ) = r.
Pour le calcul de la densité de (R, Θ), considérons une fonction h : R2 −→ R
mesurable bornée ou positive

E(h(R, Θ)) = E(h(φ(X, Y )))

Z
= h(φ(x, y))f(X,Y ) (x, y)dxdy
R 2
Z
= h(φ(x, y))fX (x)fY (y)dxdy de l'independance de X et Y
R 2
Z
1 1
= h(φ(x, y)) exp(− (x2 + y 2 ))dxdy.
R 2 2π 2

Il en découle, par changement de variable

Z ∞ Z 2π
1 1
E(h(R, Θ)) = h(r, θ) r exp(− r2 )dθdr,
0 0 2π 2

donc
1 1
f(R,Θ) (r, θ) = r exp(− r2 )1R+ (r)1[0,2π] .
2π 2
Alors on déduit que R et Θ sont indépendantes, De plus fR (r) = r exp(−r2 /2)1R+ (r),
et Θ suit la loi uniforme sur [0, 2π].
L'algorithme
p de simulation peut alors être formulé :
R ←− (−2 ln(Random))
Θ ←− 2πRandom
X ←− R cos Θ
Y ←− R sin Θ.

La formule de R dans l'algorithme précédent est justié par :

R a pour densité fR d'après la proposition 2.4, alors R2 suit une loi expo-
nentielle de paramètre 21 , en eet :

11
Soit h : R −→ R mesurable bornée ou positive, posons Z = R2

E(h(Z)) = E(h(R2 ))
Z
= h(r2 )fR (r)dr
ZR
−r 2
= rh(r2 )e 2 1R+ (r)dr,
R

et le changement de variable z = r2 , on obtient :

Z −z
e 2
E(h(Z)) = h(z) 1R+ (z)dz.
R 2

Le carré de R est simulé par la méthode d'inversion par −2 ln (Random).

Le code R suivant permet de générer par exemple 100 observations d'une
variable aléatoire gaussienne centrée réduite en utilisant la méthode de Box-
Muller.
Nous proposons l'algorithme suivant :
N = 50
X < − numeric(N )
Y < − numeric(N )
for(i in 1 : N )
{U < −runif(1)
V < −runif(1)
X[i] = sqrt(−2 ∗ log(V )) ∗ cos(2 ∗ pi ∗ U )
Y [i] = sqrt(−2 ∗ log(V )) ∗ sin(2 ∗ pi ∗ U )}
plot(X,Y)

12
Figure 3  50 observations du couple de variables indépendantes N (0, 1).

3 Calcul d'une intégrale par la méthode de Monte-

Carlo.

Les méthodes de Monte-Carlo sont particulièrement utilisées pour ap-

proximer les intégrales (en particulier, pour calculer des surfaces, des vo-
lumes, etc).
Supposons que l'on veuille estimer une quantité :
Z
I= g(x1 , x2 , ..., xd )f (x1 , x2 , ..., xd )dx1 dx2 ...dxd ,
Rd
où f est une densité.
L'idée de la méthode est de présenter I sous forme d'une espérance

I = E(g(X)),

où X variable dans Rd ayant pour densité f.

La justication théorique de la méthode est la loi forte des grands nombres
qui permet de ne faire appel qu`a une réalisation d'un échantillon, c'est à
dire à la suite (Xn (ω))n pour un seul ω .
La vitesse de convergence est bien sûr un problème crucial pour maitriser
l'erreur commise en approximant la valeur souhaitée E(X) par
Sn (ω)
X̄n (ω) = ,
n

13
que l'on peut simuler. le théorème de la limite centrale précise la vitesse de
convergence.
Théorème 3.1. : limite centrale
Soit (Xn )n≥1 une suite de variables aléatoires réelles, indépendantes et de
même loi avec E(X12 ) < ∞, alors
Sn − nE(X1 )
√
nσX1
converge en loi vers une variable N (0, 1).
Démonstration. Supposons au départ que E(X1 ) = 0 et V ar(X1 ) = 1,
Sn − nE(X1 ) Sn
√ =√ ,
nσX1 n
a comme fonction caractéristique :
Sn
it √
X 1
it √n t
ϕ√
Sn (t) = E(e n ) = (E(e ))n = (ϕX1 ( √ ))n .
n n
X1 à carré intégrable, nous pouvons donc dériver deux fois sa fonction ca-
ractéristique
0
ϕX1 (0) = 0, ϕ”X1 (0) = −1
pour tout u au voisinage de 0, un développement limité de la fonction carac-
téristique donne
u2
ϕX1 (u) = 1 − + ◦(u2 ).
2
Donc
t n
Sn (t) = lim (ϕX1 ( √ ))
lim ϕ √
n→∞ n n→∞ n
u 2 u2
= lim (1 − + ◦( ))n
n→∞ 2n n
t2
= e− 2 ,
qui est la fonction caractéristique d'une loi N (0, 1). On en déduit alors la
convergence en loi de √Sn
n
vers la loi gaussienne centrée réduite.
Si E(X1 ) 6= 0 ou V ar(X1 ) 6= 1, posons pour tout i ∈ N,
Xi − E(Xi )
Yi = ,
σXi
les Yi sont des variables centrées réduites, on applique alors la première partie
de la démonstration.

14
 Par le théorème de la limite centrale nous avons :
∀ h : R −→ R mesurable bornée avec P (N ∈ D) = 0, où N est une variable
aléatoire réelle de loi N (0, 1), et D l'ensemble des point de discontinuité de
h,
√ n Z +∞
n 1X 1 x2
lim E(h( ( g(Xk ) − E(g(X1 ))))) = h(x) √ e− 2 dx,
n→∞ σ n −∞ 2π
1
σ = σg(X1 ) .
Alors pour tous réels a et b avec a < b, nous avons aussi
n Z b
σ 1X σ 1 x2
lim P (a √ ≤ g(Xk ) − E(g(X1 )) ≤ b √ ) = √ e− 2 dx.
n→∞ n n n a 2π
1

La vitesse de convergence est en √1n , ce qui n'est pas très rapide, mais c'est
parfois la seule méthode accessible, car elle ne dépend pas de la régularité de
la fonction intégrée de plus la vitesse ne change pas si nous sommes en grande
dimension. L'erreur de la méthode de Monte-Carlo est souvent présentée soit
en donnant leur écart-type, c'est à dire √σn , soit en donnant l'intervalle de
conance à 95%. De la table de la fonction de répartition de la loi gaussienne
centrée réduite, nous avons :
n
1X σ
P (| g(Xk ) − E(g(X1 ))| ≤ 1.96 √ ) ' 0.95,
n n
1

c'est à dire que nous avons un intervalle de conance à 95%

n n
1X σ 1X σ
[ g(Xk ) − 1.96 √ , g(Xk ) + 1.96 √ ].
n n n n
1 1

Il faut se souvenir que la variance σ 2 n'est pas plus connue que la so-
lution cherchée E(g(X1 )). Nous pouvons alors faire une estimation de cette
variance.

Estimation de la variance

Soient X1 , X2 , ..., Xn un échantillon de variables aléatoires réelles, indé-

pendantes et de même loi de variance σ 2 inconnue,
n n
1 X 1X
σn2 = (g(Xi2 ) − ( g(Xk ))2 ),
n−1 n
1 1

15
et c'est un estimateur sans biais de σ 2 .
L'inégalité de Chebychev donne une estimation de la probabilité que la
moyenne empirique soit  loin  de l'espérance, mais c'est une estimation
très grossière, et peut être nettement améliorée par l'inégalité de Hoeding.
Théorème 3.2. : Hoeding
Soit (Xn )n≥1 une suite de variables aléatoires indépendantes. Supposons que
pour tout i, il existe ai , bi ∈ R tels que ai ≤ Xi ≤ bi . Alors pour tout c > 0,
n n
X X −2c2
P (| Xi − E( Xi )| ≥ c) ≤ 2 exp( Pn 2
).
1 1 1 (bi − ai )

Démonstration. Nous avons pour tout h > 0 et x ∈ R,

1[0;∞[ (x) ≤ ehx ,

remplaçons x par
Pn
1 (Xi − E(Xi )) − c
Pn
1Pn1 (Xi −E(Xi ))−c≥0 ≤ eh( 1 (Xi −E(Xi ))−c) ,

prenons l'espérance
n
X Pn
P ( (Xi − E(Xi )) ≥ c) ≤ e−hc E(eh( 1 (Xi −E(Xi )) ),
1

par indépendance des Xi il en découle

n
X n
Y
P ( (Xi − E(Xi )) ≥ c) ≤ e−hc E(eh(Xi −E(Xi )) ).
1 1

Posons alors
Yi = Xi − E(Xi ),
et montrons que
h2
(bi −ai )2
E(ehYi ) ≤ e 8 .
Pour tout i ≥ 1,
ai − E(Xi ) ≤ Yi ≤ bi − E(Xi ),
que l'on peut réécrire sous la forme

Yi = αi (ai − E(Xi )) + (1 − αi )(bi − E(Xi )),

16
avec
Yi − bi + E(Xi ) bi − Xi
αi = = .
ai − bi bi − ai
D'autre part, la fonction x 7−→ ehx est convexe, alors

ehYi ≤ αi eh(ai −E(Xi )) + (1 − αi )eh(bi −E(Xi )) ,

ce qui implique

E(ehYi ) ≤ E(αi )eh(ai −E(Xi )) + E(1 − αi )eh(bi −E(Xi )) ,

bi − E(Xi ) h(ai −E(Xi )) ai − E(Xi ) h(bi −E(Xi ))
= e − e
bi − ai bi − ai
bi − E(Xi ) hai ai − E(Xi ) hbi
= e−hE(Xi ) [ e − e ].
bi − ai bi − ai
Notons hi = h(bi − ai ),
−b
hi
bi − E(Xi ) b h−a
E(Xi ) i a ai − E(Xi ) b h−a
i b
E(ehYi ) ≤ e i −ai [ e i i i− e i i i]
bi − ai bi − ai
hi
(a −E(X )) bi − E(X i ) a i − E(X i ) hi
= e bi −ai i i
[ − e ]
bi − ai bi − ai
= eLi (hi ) ,

hi bi − E(Xi ) ai − E(Xi ) hi
avec Li (hi ) = (ai − E(Xi )) + ln( − e ).
bi − ai bi − ai bi − ai
Posons
E(Xi ) − ai
pi = ,
bi − ai
Li (hi ) = −pi hi + ln((1 − pi ) + pi ehi ).
Le développement de Taylor Lagrange d'ordre 1 de Li au voisinage de 0
donne

Li (u) = −pi u + ln((1 − pi ) + pi eu ),

0 pi eu
Li (u) = −pi + ,
1 − pi + pi eu
pi (1 − pi )eu
L”i (u) = .
((1 − pi ) + pi eu )2

17
donc
0 L”i (ξi ) 2
Li (hi ) = Li (0) + Li (0)hi + hi , ξi ∈]0, hi [
2
L”i (ξi ) 2
= hi .
2
Pour montrer
h2
(bi −ai )2
E(ehYi ) ≤ e 8 ,
il sut de vérier que
1
L”i (ξi ) ≤ .
4
pi ), en eet :
La fonction L”i possède un maximum au point ln( 1−p i

(3) pi (1 − pi )eu (1 − pi − pi eu )
Li (u) = ,
(1 − pi + pi eu )3

Li (u) = 0 donne u = ln( 1−p

pi ) qui le point critique unique de Li ,
(3) i ”

(4) pi (1 − pi )eu (1 − pi − pi eu ) − p2i (1 − pi )e2u − 3p2i (1 − pi )e2u

Li (u) =
(1 − pi + pi eu )6
pi (1 − pi )(1 − pi − 5pi eu )
= ,
(1 − pi + pi eu )3

(4) 1 − pi
Li (ln( )) = −2(1 − pi ) < 0.
pi
Ainsi
1 − pi (1 − pi )2 1
∀u ∈ R, L”i (u) ≤ L”i (ln( )) = 2
= .
pi 4(1 − pi ) 4

Nous avons donc

L”
i (ξi ) h2 1 2 h2
(bi −ai )2
E(ehYi ) ≤ eLi (hi ) = e 2 i ≤ e 8 hi = e 8 .

Il en découle alors
n n
X Y 1 2 2
−hc
P ( (Xi − E(Xi )) ≥ c) ≤ e e 8 h (bi −ai )
1 1
h2 Pn 2
= e−hc+ 8 1 (bi −ai ) .

18
Posons n
h2 X
g(h) = −hc + (bi − ai )2 ,
8
1
n
0 h X
g (h) = −c + (bi − ai )2 ,
4
1

g (h) = 0 implique
0

4c
h0 = Pn 2
1 (bi − ai )
,
n
1X
g ” (h) = (bi − ai )2 ≥ 0,
4
1

donc h0 est un minimum de g .

Finalement,
n
X h20
Pn 2
P ( (Xi − E(Xi )) ≥ c) ≤ e−h0 c+ 8 1 (bi −ai )
1
−2c2
= exp( Pn 2
).
1 (bi − ai )

Pour conclure le résultat du théorème, montrons que

n
X −2c2
P( (Xi − E(Xi )) ≤ −c) ≤ exp( Pn 2
),
1 1 (bi − ai )

pour cela, il sut de poser Zi = −Xi , puis appliquer le résultat précédent

aux variables Zi , i ∈ N.

19
4 Méthode de Monte-Carlo par Chaînes de Mar-

kov.

L'objectif est toujours de calculer I = h(x)dπ(x), où π est une loi de

R

probabilité. Dans la section précédente, on a estimé I par

n
1X
h(Xk ),
n
1

où X1 , ..., Xn est un échantillon de variables aléatoires indépendantes de

même loi π . le problème est dans la simulation d'un échantillon i.i.d de la loi
π , car elle n'est pas toujours explicitable.

La méthode de Monte-Carlo par chaîne de Markov donne une estimation

de I égale
n
1X
h(Xk ),
n
1

où X1 , ..., Xn des états successifs d'une chaîne de Markov. Avant d'aller à la

justication de la méthode, on passe par quelques rappels sur les chaînes de
Markov.

Dans tout ce chapitre, E est un ensemble dénombrable, P(E) est l'en-

semble de ses parties.

4.1 La propriété de Markov

Dénition 4.1. On dit qu'une suite de variables aléatoires (Xn )n∈N , à va-
leurs dans (E, P(E)) et dénie sur un espace probabilisé (Ω, A, P ), est une
chaîne de Markov si, pour tout (n + 1) uplet (i0 , ..., in ) de points de E tel que
P (∩n−1
j=0 {Xj = ij }) > 0,

P {Xn = in / ∩n−1
j=0 {Xj = ij }} = P {Xn = in /Xn−1 = in−1 } (1)
Autrement dit, la loi de Xn conditionnellement à (X0 , ..., Xn−1 ) et la loi de
Xn conditionnellement à Xn−1 sont identiques.
On appelle E l'espace des états. La loi de X0 est appelé la loi ou la mesure
initiale, et l'égalité (1) s'appelle propriété de Markov.

20
Dénition 4.2. On dit qu'une chaîne de Markov (Xn )n∈N est homogène si,
pour tout couple (i, j) de points de E , P (Xn+1 = j/Xn = i) est indépendant
de n, n décrivant l'ensemble des entiers pour lesquels P (Xn = i) > 0.
Dénition 4.3. Une matrice M = (Mij )i,j∈E est une matrice stochastique
si elle vérie
(i) MP ij ≥ 0 pour tout i, j ∈ E
(ii) j∈E Mij = 1 pour tout i ∈ E
Proposition 4.1. Le produit de deux matrices stochastiques est encore une
matrice stochastique.
Démonstration. Soient A = (aij )i,j∈E et B = (bij )i,j∈E deux matrices sto-
chastiques, montons que AB = (cij )i,j∈E est une matrice stochastique,
pour tout i, j ∈ E , X
cij = aik bkj ,
k∈E
on a
i- cP
ij ≥ 0,
ii- j∈E cij = j∈E k∈E aik bkj = k∈E (aik j∈E bkj ) = 1,
P P P P
donc AB est une matrice stochastique.
Soit (Xn )n∈N une chaîne de Markov homogène. On note alors pi,j la va-
leur commune des P (Xn+1 = j/Xn = i) et P = (pij )i,j∈E . La matrice P est
appelée matrice de transition de la chaîne, et c'est une matrice stochastique.

Dans la suite, X = (Xn )n≥0 est une chaîne de Markov homogène, dénie
sur un espace probabilisé (Ω, A, P), à valeurs dans (E, P(E)), de matrice de
transition P = (pij )i,j∈E et de loi initiale π0 .

4.2 Calcul des lois marginales

Théorème 4.1. :
1− Pour tout n ≥ 1 et tous i0 , ..., in ∈ E ,
P (X0 = i0 , ..., Xn = in ) = π0 (i0 )pi0 i1 ...pin−1 in .
2− Pour tous A0 , ..., An−1 , B1 , ..., Bm ∈ P(E), i ∈ E
P (Xn+1 ∈ B1 , ..., Xn+m ∈ Bm /X0 ∈ A0 , ..., Xn−1 ∈ An−1 , Xn = i)

= P (Xn+1 ∈ B1 , ..., Xn+m ∈ Bm /Xn = i)

= P (X1 ∈ B1 , ..., Xm ∈ Bm /X0 = i).

21
Démonstration . Le point 1 du théorème :
P (X0 = i0 , ..., Xn = in )

= P (Xn = in /X0 = i0 , ..., Xn−1 = in−1 )P (X0 = i0 , ..., Xn−1 = in−1 )

= P (Xn = in /Xn−1 = in−1 )P (X0 = i0 , ..., Xn−1 = in−1 ) par la propriété de Markov
= pin−1 in P (X0 = i0 , ..., Xn−1 = in−1 )
= pin−1 in P (Xn−1 = in−1 /X0 = i0 , ..., Xn−2 = in−2 )P (X0 = i0 , ..., Xn−2 = in−2 )
= pin−1 in pin−2 in−1 P (X0 = i0 , ..., Xn−2 = in−2 ),

après plusieurs itérations, nous avons

P (X0 = i0 , X1 = i1 , ..., Xn = in ) = pin−1 in ...pi0 i1 P (X0 = i0 )

= pin−1 in ...pi0 i1 π0 (i0 ).

Le point 2 :

P (Xn+1 ∈ B1 , ..., Xn+m ∈ Bm /X0 ∈ A0 , ..., Xn−1 ∈ An−1 , Xn = i)

P P P P
i0 ∈A0 ... in−1 ∈An−1 j1 ∈B1 ... jm ∈Bm pi0 i1 ...pin−1 i pij1 ...pjm−1 jm
= P P
i0 ∈A0 ... in−1 ∈An−1 pi0 i1 ...pin−1 i
X X
= ... pij1 ...pjm−1 jm
j1 ∈B1 jm ∈Bm

= P (X1 ∈ B1 , ..., Xm ∈ Bm /X0 = i).

La matrice P n = P ... × P} est une matrice stochastique, d'après la propo-

| × {z
nf ois
sition 4.1. On note pnij ses éléments.

Corollaire 4.1. Pour

P tous entiers naturels n, m et tous états i, j ∈ E,
i − P (Xn = j) = k∈E π0 (k)pnkj ,
ii − P (Xn+m = j/Xm = i) =Ppnij ,
iii− P (Xn+m = j/X0 = i) = k∈E P (Xm = j/X0 = k)P (Xn = k/X0 = i).
L'égalité (iii) est appelée équation de Chapman-Kolmogorov.

22
Démonstration. i- Elle se fait par récurrence, pour n = 1
P (X1 = j) = P (X1 = j, Ω)
X
= P (X1 = j, X0 = k)
k∈E
X
= P (X1 = j/X0 = k)P (X0 = k)
k∈E
X
= π0 (k)pkj .
k∈E

On suppose qu'elle est vraie pour n xé, et on montre qu'elle est vraie pour
n + 1,

P (Xn+1 = j) = P (Xn+1 = j, Ω)
X
= P (Xn+1 = j, Xn = i)
i∈E
X
= P (Xn+1 = j/Xn = i)P (Xn = i)
i∈E
X X
= pij π0 (k)pnki
i∈E k∈E
X
= π0 (k)pn+1
kj .
k∈E

ii- Elle se fait aussi par récurrence, pour n = 1

P (Xm+1 = j/Xm = i) = pij .

Supposons la vraie au rang n, et donc

X
P (Xm+n+1 = j/Xm = i) = P (Xm+n+1 = j, Xn+m = k/Xm = i)
k∈E
X
= P (Xm+n+1 = j/Xn+m = k, Xm = i)P (Xn+m = k/Xm = i)
k∈E
X
= pkj pnik = pn+1
ij .
k∈E

23
 X
iii − P (Xn+m = j/X0 = i) = P (Xn+m = j, Xn = k/X0 = i)
k∈E
X
= P (Xn+m = j/Xn = k)P (Xn = k/X0 = i)
k∈E
X
= P (Xm = j/X0 = k)P (Xn = k/X0 = i).
k∈E

Remarque 4.1. La donnée de la matrice de transition et de la loi initiale

sut à caractériser la loi d'une chaîne de Markov homogène jusqu'à tout
instant xé.

4.3 Récurrence et transience

Notation 4.1. :
Soit A ∈ P(E),
TA = inf{n > 0, Xn ∈ A}.
(k) (k−1)
TA = inf{n > TA , Xn ∈ A}.
.
Soient i, j ∈ E
ρij = P (Tj < ∞/X0 = i) = Pi (Tj < ∞).
X
Ni = 1{Xn =i} .
n≥1

Dénition 4.4. Soient i et j deux éléments de E. On dit que i conduit à j ,

noté i j , si ρij > 0 ; on dit que i et j communiquent, noté i ! j , si i
conduit à j et j conduit à i.

Dénition 4.5. Un point i ∈ E est dit récurrent pour la chaîne de Markov

si ρii = 1. Il est dit transient dans le cas contraire.
Dénition 4.6. Un ensemble C ⊂ E est dit clos si :
ρij = 0 ∀ i ∈ C, j ∈
/ C.

Dénition 4.7. Un ensemble C ⊂ E est dit irréductible si :

i j ouj i; ∀ i, j ∈ C.

24
Remarque 4.2. :
1− ∀ i, j ∈ E
ρij > 0 ⇐⇒ ∃n ≥ 1, pnij > 0.
En eet :
ρij > 0 ⇐⇒ Pi (Tj < ∞) > 0
⇐⇒ Pi (∪n≥1 {Xn = j}) > 0
⇐⇒ ∃n ≥ 1, Pi (Xn = j) > 0
⇐⇒ ∃n ≥ 1, pnij > 0. par le point 2 du corollaire 4.1
2− Si E est irréductible, la chaîne est dite irréductible.
3− Si tous les états sont récurrents la chaîne est dite récurrente. Elle est
dite transiente si les états sont transients.
4− La relation ” ” est une relation d'équivalence sur l'ensemble des états
récurrents, en eet :
Soient i, j, k des états récurrents. Commençons par ” ” réexive, on a
ρii = 1,

implique
i i.
Pour ” ” symétrique, montrons que si i j alors j i

i j ⇐⇒ ∃n ≥ 1, pnij > 0,

on a
pnij Pj (Ti = ∞) ≤ Pi (Ti = ∞) = 0,
car l'ensemble des trajectoires qui passent de l'états i à l'états j en n étapes,
puis qui n'atteignent jamais i est inclus dans l'ensemble des trajectoires qui
partant de i n'atteignent jamais i, ainsi
Pj (Ti = ∞) = 0,

implique
Pj (Ti < ∞) = ρji = 1,
alors j i.
Il reste à montrer ” ” est transitive, on a :
i j et j k , montrons que i k

i j ⇐⇒ ∃n ≥ 1, pnij > 0,

25
 j k ⇐⇒ ∃m ≥ 1, pm
jk > 0,

pn+m
ik ≥ pnij pm
jk > 0,

donc i → k.
Dénition 4.8. Soit i un état récurrent, on note Ei (Ti ) = E(Ti /X0 = i),
i− i est dit récurrent positif si Ei (Ti ) < ∞.
ii− i est dit récurrent nul si Ei (Ti ) = +∞.

4.4 Mesure invariante

Dénition 4.9. Soit µ une mesure positive sur E , on note aussi par µ le
vecteur ligne (µ{i}, i ∈ E). µ est dit mesure invariante de la chaîne si
µP = µ.

On prendra garde au fait que µ n'est pas nécessairement une probabilité.

Observons que si π0 est une mesure invariante de la chaîne, Xn est de loi π0
pour tout n ∈ N.
Notation 4.2. On note ER l'ensemble des états récurrent, ET l'ensemble
des états transitoires.
n
X
Nn (j) = 1{Xi =j} , j ∈ E
i=1

mij = Ei (Tj ), i, j ∈ E

Théorème 4.2. Soit j ∈ ER , on a

Nn (j) 1
lim = 1 p.s.
n→∞ n mjj {Tj <∞}

Pour i ∈ E , on a
Nn (j) ρij
lim Ei ( )= .
n→∞ n mjj
Théorème 4.3. Soit (Xn )n≥0 une chaîne de Markov homogène, irréductible
récurrente positive de matrice de transition P , alors (Xn )n≥0 admet une loi
invariante π, donnée par :
1
π({i}) = .
mii

26
Démonstration. Supposons d'abord que (Xn )n≥0 admet une loi invariante π,
et montrons que
1
π(i) = , ∀i ∈ E
mii
π loi invariante, elle vérie :
π = πP,
c'est-à-dire pour tout i ∈ E ,
X
π(i) = π(j)pji ,
j∈E
ainsi X XX
1= π(i) = π(j)pji ,
i∈E i∈E j∈E
on a
n n
X 1 X m 1 XX
π(j)( pji ) = π(j)pm
ji
n n
j∈E m=1 m=1 j∈E
n X
1 X
= π(j)P (Xm = i/X0 = j)
n
m=1 j∈E
n X
1 X
= P (Xm = i, X0 = j)
n
m=1 j∈E
n
1 X
= π(i)
n
m=1
= π(i),
or
n n
1 X m 1 X
lim pji = lim P (Xm = i/X0 = j)
n→∞ n n→∞ n
m=1 m=1
n
1 X
= lim Pj (Xm = i)
n→∞ n
m=1
n
1 X
= lim Ej (1Xm =i )
n→∞ n
m=1
Nn (i)
= lim Ej ( )
n→∞ n
ρji
= , d'après le théorème 4.2.
mii

27
ρji = 1, en eet :
(Xn )n≥0 chaîne de Markov irréductible, alors

∀i, j ∈ E, i j ou j i,

puisque la chaîne est récurrente, et ” ” symétrique sur l'ensemble des

états récurrents, on a

∀i, j ∈ E, i j et j i,

i j ⇐⇒ ∃n ≥ 1, pnij > 0,

on a
pnij Pj (Ti = ∞) ≤ Pi (Ti = ∞) = 0,
ainsi
Pj (Ti = ∞) = 0,
implique
Pj (Ti < ∞) = ρji = 1
.
Par passage à la limite quand n tend vers l'inni dans l'équation
n
X 1 X m
π(j)( pji ) = π(i),
n
j∈E m=1

on trouve X 1
π(j) = π(i),
mii
j∈E

implique
1
π(i) = .
mii
Montrons maintenant l'existence de la loi invariante, soit E1 ⊂ E ni, on a :
X 1 Xn n
1 X X m
( pm
ji ) = ( pji )
n n
i∈E1 m=1 m=1 i∈E1
n X
1 X
≤ ( pm
ji )
n
m=1 i∈E
= 1,

28
on prend la limite quand n tend vers l'inni, on obtient
X 1
≤ 1,
mii
i∈E1

comme E1 est arbitraire, on a

X 1
≤ 1. (2)
mii
i∈E

D'autre part,
X 1 Xn n
m 1 X X m
( pki )pij = ( pki pij )
n n
i∈E1 m=1 m=1 i∈E1
n X
1 X
≤ ( pm
ki pij )
n
m=1 i∈E
n
1 X
= pm+1
kj ,
n
m=1

quand n −→ ∞, on a
X 1 1
pij ≤ ,
mii mjj
i∈E1

comme E1 est arbitraire

X 1 1
pij ≤ , (3)
mii mjj
i∈E

en sommant sur j ,
X 1 X 1
≤ ,
mii mjj
i∈E j∈E
or X 1 X 1
= ,
mii mjj
i∈E j∈E

donc les inégalités précédentes sont des égalités, on pose π(i) = mii ,
1

X
π(i)pij = π(j),
i∈E

d'où l'existence de la loi invariante.

29
Dénition 4.10. Soit π une loi de probabilité sur E . La matrice de transi-
tion P est dite réversible par rapport à π si
π({i})pij = π({j})pji ,

pour tout i, j ∈ E .
Proposition 4.2. Si P est réversible par rapport à π alors π est une proba-
bilité invariante.
Démonstration. Pour montrer que π est une probabilité invariante, il sut
de montrer que : X
∀i ∈ E, π({i}) = π({j})pji .
j∈E

π({i})pij car π est réversible

X X
π({j})pji =
j∈E j∈E
X
= π({i})P (X1 = j/X0 = i)
j∈E

= π({i}).

4.5 Période d'une chaîne de Markov

Dénition 4.11. Soit i ∈ E tel que ρii > 0, on pose

di = pgcd{n ≥ 1/ pnii > 0}

di est appelé période de i.

Théorème 4.4. Soient i, j ∈ E
Si i ! j alors di = dj .
Démonstration. Soient i, j ∈ E tel que i ! j , on a ρij > 0 et ρji > 0, alors
∃k ≥ 1, pkij > 0,

et
∃l ≥ 1, plji > 0,
ainsi
pk+l k l
ii ≥ pij pji > 0,

30
alors di divise k + l.
Soit n ≥ 1 tel que pnjj > 0
pk+l+n
ii ≥ pkij pnjj plji > 0,

donc di divise k + l + n.
On a di divise k + l et divise k + l + n, donc divise n, ce qui implique que
di ≤ dj ,

de même on montre que

dj ≤ di ,
d'où
di = dj .

Conséquence 4.1. Si (Xn )n≥0 est une chaîne de Markov irréductible récur-
rente, alors
∀i, j ∈ E di = dj = d.
d est appelé la période de la chaîne, si d = 1 la chaîne est dite apériodique.

4.6 Méthode de Monte-Carlo

Soit E un espace d'état ni et π la loi invariante de matrice de tran-

Rsition P irréductible et apériodique. Supposons que l'on cherche à calculer
h(x)dπ(x), où h est une fonction donnée. Une méthode naturelle dite de
Monte Carlo par chaîne de Markov consiste à générer les états successifs
X0 ...Xn d'une chaîne de Markov de matrice de transition P à partir d'un
état initial X0 de loi quelconque. Dans le théorème suivant, on établit que
n
1X
Sn (h) = h(Xk )
n
1

est une bonne approximation de h(x)dπ(x), et la loi de Xn est proche de

R

π.
Théorème 4.5. Soit (Xn )n≥0 une chaîne de Markov sur E irréductible ré-
currente positive, π la loi invariante
R de la chaîne alors pour toute fonction
f dénie sur E telle que f ≥ 0 ou f (x)dπ(x) < ∞, on a :
n Z
1X
f (Xk ) −→ f (x)dπ(x) p.s
n
1

31
Algorithme de Metropolis

Dans de nombreux cas, la loi π que l'on cherche à simuler n'est pas donnée
a priori comme la mesure invariante d'une chaîne de Markov. l'algorithme
de Metropolis produit une chaîne de Markov réversible par rapport à π .
On se donne une matrice de transition markovienne Q sur E , appellée matrice
de sélection, telle que pour tout couple (i, j) ∈ E

Qij > 0 =⇒ Qji > 0

Soit h :]0, ∞[−→]0, 1[ une fonction vériant

1
h(u) = uh( ).
u
Par exemple
h(u) = inf(1, u),
ou bien
u
h(u) = .
1+u
Pour i 6= j , posons
(
Rij = π(i)Q(i,j) ) si Q(i, j) 6= 0
h( π(j)Q(j,i)
0 sinon

On construit alors une matrice de transition P dénie par

pij = Qij Rij

pour i 6= j , et X
pii = 1 − pij
i6=j

Proposition 4.3. On suppose que π charge tous les points de E . La matrice

P dénie précédemment est réversible par rapport à π . Elle est irréductible
si Q est irréductible. Si de plus h(u) < 1, elle est apériodique.
L'algorithme swivant, correspondant à la simulation de la chaîne de ma-
trice de transition P , s'appelle l'algorithme de Metropolis.
Etape 0.Initialiser X0 .
Etape n+1.
(sélection) choisir y avec la loi Q(Xn , y).
Tirer un nombre U au hasard dans [0, 1].

32
Si U < R(Xn , y) accepter la sélection.
Xn+1 = y.
Sinon refuser la sélection
Xn+1 = Xn .

5 Conclusion

Comme conclusion à cette modeste étude, nous retenons que la méthode

de Monte-Carlo est une méthode d'approximation par introduction de procé-
dés aléatoires. Cela permet d'estimer des valeurs numériques et de caractéri-
ser des systèmes complexes. L'inconvénient est que cette méthode est lente,
mais il existe des cas où c'est la seul technique accessible, en eet elle ne
dépend pas de la régularité de la fonction à intégrer. D'autre part ne dé-
pend pas de la dimension de l'espace en question ce qui évite les problèmes
de déreglement de la convergence de l'estimateur lorsqu'il sagit de grandes
dimensions.

33
6 Annexe

Théorème 6.1. :loi forte des grands nombres

Soit (Xn )n≥1 une suite de variables aléatoires indépendantes, identiquement
distribuées intégrables. Alors Snn converge presque sûrement vers E[X1 ], avec
Sn = X1 + ... + Xn .

Démonstration. La démonstration consiste à prouver dans un premier temps

le résultat sous l'hypothèse plus forte E(|X1 |4 ) < ∞ et E(X1 ) = 0. Dans ce
cas, dont on peut se contenter en première lecture, P (|Sn |/n ≥ ) peut être
majoré en utilisant l'inégalité de Markov. La borne ainsi obtenue est le terme
général d'une série convergente, ce qui permet de conclure grâce au lemme
de Borel-Cantelli. Sous l'hypothèse plus faible du théorème, on approxime
toute variable de L1 par des variables de L4 , puis on se ramène au cas
traité. Commençons donc par montrer le résultat lorsque E(|X1 |4 ) < ∞ et
E(X1 ) = 0. Dans ce cas, l'inégalité de Markov montre que pour tout n ≥ 1
et tout δ > 0,
1
P (|Sn | ≥ δn) ≤ E(Sn4 ).
δ 4 n4
Observons que
Xn
Sn4 = ( Xi )4
1
Xn X
=( Xi2 + Xi Xj )2
1 i6=j
Xn n
X X X
=( Xi2 )2 + 2 Xi2 Xj Xk + ( Xi Xj )2
1 1 j6=k i6=j
n
X X X X
= Xi4 + Xi2 Xj2 + 2( Xi Xj Xk2 + 2 Xi3 Xj )
1 i6=j 1≤i,j,k distincts≤n 1≤i6=j≤n
X X X
+2 Xi2 Xj2 +4 Xi Xj Xk2 Xi Xj Xk Xl
i6=j 1≤i,j,k distincts≤n 1≤i,j,k,l distincts≤n
X X X X
= Xi4 +4 Xi3 Xj +3 Xi2 Xj2 + 6 Xi Xj Xk2
1≤i≤n 1≤i6=j≤n 1≤i6=j≤n 1≤i,j,k distincts≤n
X
+ Xi Xj Xk Xl .
1≤i,j,k,l distincts≤n

34
Donc, par linéarité de l'espérance, indépendance et centrage des Xi ,
X X X
E(Sn4 ) = E(Xi4 ) + 4 E(Xi3 )E(Xj ) + 3 E(Xi2 )E(Xj2 )
1≤i≤n 1≤i6=j≤n 1≤i6=j≤n
X X
+6 E(Xi )E(Xj )E(Xk2 ) + E(Xi )E(Xj )E(Xk )E(Xl )
1≤i,j,kdistincts≤n 1≤i,j,k,ldistincts≤n

= nE(X14 ) + 3n(n − 1)(E(X12 ))2 .

La série
X [nE(X 4 ) + 3n(n − 1)(E(X 2 ))2 ] E(X14 ) X 1 (E(X12 ))2 X n − 1
1 1
= +3
n4 δ 4 δ4 n3 δ4 n3
n≥1 n≥1 n≥1

est convergente, ce qui donne la convergence de la serie n≥1 P (|Sn | > nδ).
P

D'aprés la proposition suivante nous avons la convergence de Sn

n vers 0
presque sûrement.
Proposition 6.1. Soient (Xn )n∈N , X des variables aléatoires réelles.
Si ∀ > 0 X
P (|Xn − X| > ) < ∞
n≥0
alors Xn −→ X p.s.
Démonstration. Posons :
An = {|Xn − X| > }.
D'après Borel-Cantelli
P (lim sup An ) = 0,
n→∞
i.e
P (∀n ≥ 0, ∃k ≥ n/|Xn − X| < ) = 1,
ce qui implique que Xn −→ X p.s .
Supposons maintenant que les Xi sont intégrables et centrées, sans autre
hypothèses. Soit  > 0 xé. Il existe, pour tout i ≥ 1, des variables Yi étagées,
centrées,P
indépendantes et de même loi, telles que E(|Xi − Yi |) ≤ .
Si Tn = 1≤i≤n Yi , nous avons
n n
1 1 X X
|Sn | = | (Xi − Yi ) + Yi |
n n
1 1
1 X 1
≤ |Xi − Yi | + |Tn |.
n n
1≤i≤n

35
Par le point précédent Tnn converge p.s. vers 0, car les Yi sont indépendantes,
de même loi , centrées et étagées ( étagées donc dansL4 ). Donc il sut de
montrer que :
1 X
lim sup |Xi − Yi |,
n→∞ n
1≤i≤n

peut être rendu arbitrairement petit en prenant  arbitrairement petit.

Posons Zi = |Xi − Yi |, i ≥ 1, soient k ∈ N et δ > 0,
1 X
P( max Zi ≥ 2E(Z1 ) + δ)
2k <n≤2k+1 n
1≤i≤n

1 X
= P( max Zi ≥ 2E(Z1 ) + δ, (∪1≤i≤2k+1 {Zi ≤ 2k }) ∪ (∪1≤i≤2k+1 {Zi > 2k }))
2k <n≤2k+1 n
1≤i≤n
1 X
≤ P ( max Zi ≥ 2E(Z1 ) + δ, ∪1≤i≤2k+1 {Zi ≤ 2k })
2k <n≤2k+1 n
1≤i≤n
1 X
+ P ( max Zi ≥ 2E(Z1 ) + δ, ∪1≤i≤2k+1 {Zi > 2k })
2k <n≤2k+1 n
1≤i≤n
1 X
≤ P ( max Zi ≥ 2E(Z1 ) + δ, ∪1≤i≤2k+1 {Zi ≤ 2k }) + P (∪1≤i≤2k+1 {Zi > 2k })
2k <n≤2k+1 n
1≤i≤n
1 X
= P ( max Zi 1[0,2k ] (Zi ) ≥ 2E(Z1 ) + δ) + P (∃i ∈ {1, 2, ..., 2k+1 } : Zi > 2k ).
k
2 <n≤2 k+1 n
1≤i≤n

Nous avons
1 X 1 X
max Zi 1[0,2k ] (Zi ) ≤ max max Zi 1[0,2k ] (Zi )
2k <n≤2k+1 n 2k <n≤2k+1 n 2k <n≤2k+1
1≤i≤n 1≤i≤n
1 X
= k Zi 1[0,2k ] (Zi ),
2 k+1
1≤i≤2

Ce qui implique
1 X 1 X
{ max Zi 1[0,2k ] (Zi ) ≥ 2E(Z1 )+δ} ⊂ { k Zi 1[0,2k ] (Zi )) ≥ 2E(Z1 )+δ}.
2k <n≤2k+1 n 2
1≤i≤n 1≤i≤2k+1

De plus on a
k+1
2X k+1
2X
k+1 k
{ Zi 1[0,2k ] (Zi ) ≥ 2 E(Z1 )+2 δ} ⊂ { Zi 1[0,2k ] (Zi ) ≥ 2k+1 E(Z1 1[0,2k ] (Z1 ))+2k δ}
1 1

36
Donc
2 k+1
1 X
P( max Zi > 2E(Z1 ) + δ)
2k <n≤2k+1 n
1

X
≤ 2k+1 P (Z1 > 2k ) + P ( Zi 1[0,2k ] (Zi ) ≥ 2k+1 E(Z1 ) + 2k δ)
1≤i≤2k+1
X
≤ 2k+1 P (Z1 > 2k ) + P ( [Zi 1[0,2k ] (Zi ) − E(Zi 1[0,2k ] (Zi ))] ≥ 2k δ).
1≤i≤2k+1

Par l'inégalité de Tchebitchev

X
P( [Zi 1[0,2k ] (Zi ) − E(Zi 1[0,2k ] (Zi ))] ≥ δ2k )
1≤i≤2k+1
− E(Zi 1[0,2k ] (Zi ))])2
P
E( 1≤i≤2k+1 [Zi 1[0,2k ] (Zi )
≤
22k δ 2
2
= [E(Z1 1[0,2k ] (Z1 ))2 − E 2 (Z1 1[0,2k ] (Z1 )]
2k δ 2
2
≤ E(Z12 1[0,2k ] (Z1 )).
2k δ 2
Alors
1 X 2
P( max Zi ≥ 2E(Z1 )+δ) ≤ 2k+1 P (Z1 > 2k )+ k 2 E(Z12 1[0,2k ] (Z1 )).
2k <n≤2k+1 n 2 δ
1≤i≤n

Pour t < 2k+1 ,

P (Z1 > t) ≥ P (Z1 > 2k+1 ),
implique
Z 2k+1
P (Z1 > t)dt ≥ 2k P (Z1 > 2k+1 ),
2k
ce qui donne
X X
2k P (Z1 > 2k+1 ) = 2k−1 P (Z1 > 2k )
k≥−1 k≥0
≤ E(Z1 ),

multiplions les deux cotés par 4

X
2k+1 P (Z1 > 2k ) ≤ 4E(Z1 ).
k≥0

37
De plus
X X
2−k E(Z12 1[0,2k ] (Z1 )) = E(Z12 2−k 1[0,2k ] (Z1 ))
k≥0 k≥0
≤ 4E(Z1 ).

En eet :
pour chaque ω soit 0 ≤ Z1 (ω) ≤ 1 ou bien ∃p ≥ 0/ 2p < Z1 (ω) ≤ 2p+1 ,
pour 2p < Z1 (ω) ≤ 2p+1
X X
Z12 (ω) 2−k 1[0,2k ] (Z1 (ω)) = Z12 (ω) 2−k 1[0,2k ] (Z1 (ω))
k≥0 k≥p+1
X
≤ 22p+2 2−k = 2p+2
k≥p+1

≤ 4Z1 (ω),

de même pour 0 ≤ Z1 (ω) ≤ 1, donc ∀ω

X
Z12 (ω) 2−k 1[0,2k ] (Z1 (ω)) ≤ 4Z1 (ω),
k≥0

Ce qui implique le résultat cherché.

Finalement,
X 1 X
P( max Zi ≥ 2E(Z1 ) + δ) ≤ 4(1 + 2δ −2 )E(Z1 )
2k <n≤2k+1 n
k≥0 1≤i≤n

< ∞,

d'aprés le lemme de Borel-Cantelli

1 X
P (∃k ≥ 0, ∀p ≥ k/ max Zi ≥ 2E(Z1 ) + δ) = 0,
2p <n≤2p+1 n
1≤i≤n

alors pour k assez grand, nous avons :

1 X
max Zi < 2E(Z1 ) + δ p.s,
2k <n≤2k+1 n
1≤i≤n

implique pour n assez grand,

1 X
Zi < 2E(Z1 ) + δ p.s,
n
1≤i≤n

38
puisque δ est arbitraire,
1 X
lim sup Zi ≤ 2E(Z1 ).
n→∞ n
1≤i≤n

Nous pouvons maintenant nir la démonstration

1 1 X 1
lim sup |Sn | ≤ lim sup |Xi − Yi | + lim sup |Tn |
n→∞ n n→∞ n n→∞ n
1≤i≤n
≤ 2E(|X1 − Y1 |
≤ 2.

Puisque  est arbitraire, ceci conclut la démonstration.

39
Références

[1] M. Benaim, N. El Karaui. Promenade aléatoire. Chaines de

Markov et simulation, martingales et stratégies. Ecole Poly-
technique. Novembre 2006.
[2] P. Barbe, M. Ledoux. Probabilité. EDP Sciences. 2007.
[3] C. Edouard Pster. Théorie des probabilités, cours d'introduc-
tion avec application à la statistique mathématique. Presses
polytechniques et universitaires romandes. 2012
[4] L. Elie, B. Lapeyre. Introduction aux Méthodes de Monte-
Carlo. Lien : cermics.enpc.fr/ bl/PS/SIMULATION-X/poly-
monte-carlo-x.pdf. Septembre 2001.
[5] J. François Delmas, B. Jourdain. Modèles aléatoires. Appliation
aux sciences de l'ingénieur et du vivant. Springer. 2006.
[6] M. Ledra. La méthode Monté Carlo et ses application.
Lien : https ://www.researchgate.net/publication/310828308.
Mai 2016
[7] T. Mourid. Cours de M2 probabilités statistiques. Module
"Processus stochastique". 2016/2017.
[8] B. Ycart. Modèles et Algorithmes Markoviens. Springer. 2000.

40