Année 2021-2022

Simulation Monte Carlo

Madalina Deaconu

Institut Élie Cartan & Inria Nancy - Grand Est

Université de Lorraine
Email : [Link]@[Link]
Table des matières

1 Introduction 2

2 Introduction 3

3 Simulation de nombres aléatoires 4

3.1 Simulation d’une loi uniforme sur [0, 1] . . . . . . . . . . . . . . . . . . . . . . 4

4 Simulation de variables aléatoires 6

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
4.2 Simulation de lois classiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
4.2.1 Loi de Bernoulli de paramètre p . . . . . . . . . . . . . . . . . . . . . . 7
4.2.2 Loi binomiale de paramètres (n, p) . . . . . . . . . . . . . . . . . . . . 8
4.2.3 Loi uniforme sur {0, 1, ..., n − 1} . . . . . . . . . . . . . . . . . . . . . . 8
4.2.4 Loi uniforme sur [a, b] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.2.5 Loi exponentielle de paramètre λ . . . . . . . . . . . . . . . . . . . . . 10
4.2.6 Loi géométrique de paramètre p . . . . . . . . . . . . . . . . . . . . . . 10
4.2.7 Loi de Poisson de paramètre λ . . . . . . . . . . . . . . . . . . . . . . . 11
4.2.8 Loi gaussienne centrée réduite : Méthode de Box-Muller . . . . . . . . . 13
4.3 Méthodes générales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.3.1 Méthode générale pour une variable aléatoire discrète . . . . . . . . . . 14
4.3.2 Méthode de simulation par inversion de la fonction de répartition . . . 15
4.3.3 Algorithme par rejet . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.3.4 Simulation par composition . . . . . . . . . . . . . . . . . . . . . . . . 18
4.4 Simulation de vecteurs aléatoires . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.4.1 Cas indépendant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.4.2 Vecteur gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.4.3 Cas général . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

1
Chapitre 1

Introduction

L’objectif de ce cours est d’introduire les outils théoriques et d’illustrer, à travers des
problèmes réels, comment la modélisation stochastique et les méthodes de Monte Carlo per-
mettent d’apporter des réponses pour des questions issues de nombreux domaines applicatifs.
La modélisation probabiliste est fondamentale dans tous les domaines d’application.
Les outils probabilistes que nous développons, comme les méthodes de Monte Carlo et le
mouvement brownien, sont des outils généraux utilisés dans de nombreux domaines comme en
physique (écoulement d’un fluide, trajectoire d’un avion), en biologie (mutation du genôme),
en chimie (coagulation des polymères), en climatologie (modélisation du vent, de la pluie
etc), en informatique, en médicine, en finance (évaluation des produits dérivés), en science
du vivant, en assurance, en linguistique, en sociologie... Les modèles que nous étudions sont
inspirées du monde de la finance, des files d’attente, etc.
Nous débuterons ce cours par des rappels sur les méthodes numériques permettant de
simuler des variables aléatoires.
Ensuite nous étudions des méthodes numériques basées sur les tirages successifs de nombres
aléatoires, méthodes appelées généralement méthodes de Monte Carlo.
En outre, nous aborderons l’étude de certains processus stochastiques essentiels dans la
modélisation, tels que le mouvement brownien. Ce processus est fondamental en finance.

2
Chapitre 2

Introduction

L’objectif de ce cours est d’introduire les outils théoriques et d’illustrer, à travers des
problèmes réels, comment la modélisation stochastique et les méthodes de Monte Carlo per-
mettent d’apporter des réponses pour des questions issues de nombreux domaines applicatifs.
La modélisation probabiliste est fondamentale dans tous les domaines d’application.
Les outils probabilistes que nous développons, comme les méthodes de Monte Carlo et le
mouvement brownien, sont des outils généraux utilisés dans de nombreux domaines comme en
physique (écoulement d’un fluide, trajectoire d’un avion), en biologie (mutation du genôme),
en chimie (coagulation des polymères), en climatologie (modélisation du vent, de la pluie
etc), en informatique, en médicine, en finance (évaluation des produits dérivés), en science
du vivant, en assurance, en linguistique, en sociologie... Les modèles que nous étudions sont
inspirées du monde de la finance, des files d’attente, etc.
Nous débuterons ce cours par des rappels sur les méthodes numériques permettant de
simuler des variables aléatoires.
Ensuite nous étudions des méthodes numériques basées sur les tirages successifs de nombres
aléatoires, méthodes appelées généralement méthodes de Monte Carlo.
En outre, nous aborderons l’étude de certains processus stochastiques essentiels dans la
modélisation, tels que le mouvement brownien. Ce processus est fondamental en finance.

3
Chapitre 3

Simulation de nombres aléatoires

Nous allons voir dans ce qui suit que, à la base de toute méthode de simulation d’une
variable aléatoire, se trouve la loi uniforme. En outre, nous montrerons comment on peut
construire des générateurs de nombres aléatoires.

3.1 Simulation d’une loi uniforme sur [0, 1]

Dans tous les logiciels classiques dédiés au calcul numérique, il existe de très bons générateurs
de nombres aléatoires.
Nous présentons leur principe général de fonctionnement.
La méthode la plus couramment utilisée est la méthode des congruences linéaires.
On simule un nombre aléatoire sur [0, 1] comme suit :
On construit une suite (xn )n≥0 de nombres entiers compris entre 0 et m − 1 de la façon
suivante :
1. On choisit un nombre entier quelconque, appelé valeur initiale

x0 = valeur initiale. (3.1)

2. On construit la suite xn par récurrence

xn+1 = axn + b (modulo m), (3.2)

avec a, b, m des entiers positifs, qu’il faudra choisir soigneusement, si l’on veut que les
caractéristiques de la suite soient performantes.
Sedgewick préconise le choix suivant :

 a = 31415821
b=1
m = 108 .


Cette méthode permet de simuler des entiers pseudo aléatoires entre 0 et m − 1.

3. Pour obtenir un nombre aléatoire entre 0 et 1 on divise le résultat, cet entier aléatoire
ainsi généré, par m
xn
u= . (3.3)
m

4
On dit que de tels nombres sont (pseudo-)aléatoires pusiqu’ils sont obtenus par une règle
arithmétique. Mais la grande périodicité de m et un choix judicieux de a et b permettent
d’affirmer (cf. Knuth) que leur loi de distribution est ”indistinguable presque sûrement” de
vrais nombres aléatoires.
Le générateur précédent fournit des résultats acceptables dans les cas courants. Cepen-
dant, sa période (ici m = 108 ) peut se réveler parfois insuffisante.
On peut alors, obtenir des générateurs de nombres aléatoires de période arbitrairement
longue en augmentant m.

5
Chapitre 4

Simulation de variables aléatoires

4.1 Introduction
On présentera dans cette partie un panel de méthodes pour la simulation de variables
aléatoires. Nous commençons par rappeler les principes de simulation pour les lois classiques.
Ensuite nous introduisons les méthodes générales de simulation pour les variables aléatoires
à valeurs discrètes et continues.
Soit (Ω, F, P) un espace de probabilité et X une variable aléatoire réelle sur cet espace.
On notera F sa fonction de répartition et f sa densité. Plus précisément, pour tout x ∈ R :

F (x) = P(X ≤ x)

et
d
F (x).
f (x) =
dx
Principe Pour simuler des variables aléatoires de loi quelconque, l’idée est de se ramener
à la loi uniforme sur [0, 1]. La loi uniforme sur [0, 1] est donc à la base de toute simulation
de variable aléatoire. Rappelons ses propriétés.

Soit U une variable aléatoire de loi uniforme sur [0, 1], U ∼ U[0, 1]. Dans ce cas :

 0 si x < 0
FU (x) = x si 0 ≤ x ≤ 1
1 si x > 1


et 
1 si x ∈ [0, 1]
fU (x) = 1{x∈[0,1]} =
0 si x ∈
/ [0, 1].
Cadre général : On suppose qu’on dispose d’un générateur de variables aléatoires de loi
uniforme sur [0, 1] indépendantes.
L’hypothèse d’indépendance des valeurs est une des conditions essentielles pour la validité
de la plupart des algorithmes présentés dans la suite.

Exemple 1. Matlab possède une fonction rand() dont les appels successifs fournissent une
suite de variables aléatoires “indépendantes” et identiquement distribuées, de loi uniforme
sur [0, 1]. Nous ne nous intéresserons pas ici à la conception d’une telle fonction.

6
 Les générateurs sont souvent construits à partir d’une relation de congruence sur de
nombres de grande dimension et initialisés par exemple à partir de l’horloge de la machine. Il
s’agit des générateurs de nombres pseudo-aléatoires, notamment les valeurs obtenues ne sont
qu’apparemment indépendantes.
Remarque 1. Une variable aléatoire U de loi uniforme sur [0, 1] est différente de 0 et 1
presque surement. En particulier, si Λ(dx) note la mesure de Lebesgue, on a
Z
P(U ∈ {0, 1}) = 1[0,1] dΛ = Λ({0, 1}) = 0.
{0,1}

Remarque 2. En Matlab la fonction rand ne renvoie que des valeurs différentes de 0 et 1.

Objectif du cours : Vous apprendre à construire des méthodes pour simuler une variable
aléatoire ou un vecteur aléatoire suivant une loi donnée, différents de la loi uniforme sur [0, 1].

4.2 Simulation de lois classiques

4.2.1 Loi de Bernoulli de paramètre p
Soit p ∈ [0, 1]. On veut simuler une variable aléatoire X ∼ B(p) de loi de probabilité
µ(x) = pδ1 (x) + (1 − p)δ0 (x). (4.1)
On peut exprimer cela aussi sous la forme X est une variable aléatoire à valeurs dans {0, 1}
et sa loi est donnée par p0 = P(X = 1) = p et p1 = P(X = 0) = 1 − p. On remarque que
p0 + p1 = 1.

Algorithme Pour simuler X on tire U une uniforme sur [0, 1] et on définit


1 si U ≤ p
X=
0 sinon,
c’est-à-dire X = 1{U ≤p} .
On propose deux fonctions pour simuler cette variable aléatoire. La première ne rend
qu’une réalisation d’une variable aléatoire de la loi de Bernoulli de paramètre p, la seconde
renvoie un k-échantillon, qui sera utilisé par la suite dans la simulation de la loi binomiale :

Algorithmes Matlab
function x=bernoulli1(p)
% tirage suivant la loi de Bernoulli de parametre p
x=0;
if rand()<p
x=1;
end

function x=bernoulli2(k,p)
% tirage d’un k-echantillon suivant la loi de Bernoulli de parametre p
y=rand(1,k),
x=(y<p);
end

7
Exemple 2. Pile ou face
On veut simuler une variable aléatoire de loi “Pile ou face” i.e. X ∼ B( 12 ). On tire U une
uniforme sur [0, 1] et on définit
( 1
“Pile” si U ≤
X= 2
“Face” sinon,

formellement on peut aussi écrire X = 1{U ≤ 1 } .

2

4.2.2 Loi binomiale de paramètres (n, p)

Soient n ∈ N∗ et p ∈ [0, 1]. On veut simuler X ∼ B(n, p) donc de loi
n  
X n k
µ(x) = p (1 − p)n−k δk (x). (4.2)
k=0
k

Donc X est une variable aléatoire à valeurs dans {0, . . . , n} de loi donnée par les probabilités
pk = P(X = k) = nk pk (1 − p)n−k ,pour tout k ∈ {0, . . . , n}.

On sait qu’une variable aléatoire binomiale peut être représentée comme la somme de n
variables aléatoires indépendantes de Bernoulli de paramètre p. Plus précisément on a le
résultat suivant :

Lemme 1. Soient (Yi )1≤i≤n des variables aléatoires indépendantes

Pn et identiquement dis-
tribuées de loi de Bernoulli de paramètre p, alors X = i=1 Yi suit une loi binomiale de
paramètres (n, p).

Algorithme
Il suffit donc de simuler n variables aléatoires indépendantes de loi B(p) et d’en faire la
somme. On simule (Ui )1≤i≤n , n v.a.i.i.d. de loi uniforme sur [0, 1] et on définit
n
1{Ui ≤p} .
X
X=
i=1

Algorithme Matlab
Nous utilisons ici la fonction Matlab définie précédemment qui simule un n échantillon de
loi de Bernoulli de paramètre p.

function x=binomiale(n,p)
% tirage suivant la loi binomiale de parametres (n,p)
x=sum(bernoulli2(n,p));

4.2.3 Loi uniforme sur {0, 1, ..., n − 1}

On veut simuler X de loi uniforme sur {0, 1, ..., n − 1}, donc de loi
n−1
1X
µ(x) = δk (x). (4.3)
n k=0

Méthode 1 : Nous allons utiliser le résultat suivant :

8
Lemme 2. Si U est une variable aléatoire de loi uniforme sur [0, 1] alors la partie entière
X = [nU ] suit la loi uniforme sur {0, 1, ..., n − 1}.
Remarque 3. En admettant ce résultat, l’algorithme s’écrit :
Algorithme Matlab
function x=uniforme1(n)
% tirage suivant la loi uniforme sur {0,...,n-1}
x=floor(n*rand());
end
Démonstration : On note X la variable aléatoire rendue par cette fonction.
Comme P(0 ≤ rand() < 1) = 1, on a P(0 ≤ n*rand() < n) = 1 et
P(floor(n*rand()) ∈ {0, 1, ..., n − 1}) = 1.
Donc X prend ses valeurs dans {0, 1, ..., n − 1}.
Maintenant, soit k ∈ {0, 1, ..., n − 1} :
P(X = k) = P(floor(n*rand()) = k) = P(k ≤ n*rand() < k + 1)
 
k k+1 1
= P ≤ rand() < = .
n n n
Ce qui prouve le résultat souhaité. 

Méthode 2 : En utilisant la méthode générale de simulation d’une variable aléatoire discrète,

voir 4.3.1.

4.2.4 Loi uniforme sur [a, b]

On souhaite simuler X une variable aléatoire de loi uniforme sur l’intervalle [a, b] pour
a < b, avec densité
( 1
1
f (x) = 1[a,b] (x) = b − a si x ∈ [a, b] (4.4)
b−a 0 si x ∈
/ [a, b].
Lemme 3. Si U suit la loi uniforme sur [0, 1], alors, si a < b et a, b ∈ R, la variable aléatoire
X = a + (b − a)U suit la loi uniforme sur [a, b].
Démonstration : Soit ϕ : R → R une fonction continue bornée. En faisant le changement
de variable x = a + (b − a)u, on a
Z 1 Z b
1
E(ϕ(a + (b − a)U )) = ϕ(a + (b − a)u)du = ϕ(x) dx,
0 a b−a
ce qui signifie que X est une variable aléatoire de densité f (x) = 1
b−a
1[a,b] (x), et on reconnaı̂t
ainsi la densité de la loi uniforme sur [a, b]. 
L’algorithme associé s’écrit alors :

Algorithme Matlab :
function x=uniforme2(a,b)
% tirage suivant la loi uniforme sur [a,b]
x=a+(b-a)*rand();
end

9
4.2.5 Loi exponentielle de paramètre λ
Soit λ > 0. On souhaite simuler une variable aléatoire X de loi exponentielle de paramètre
λ, avec 
λ exp(−λx) si x ≥ 0
f (x) = λ exp(−λx)1R+ (x) = (4.5)
0 si x < 0.
Sa fonction de répartition est F (x) = P(X ≤ x) = 1 − exp(−λx). Cette fonction est une
bijection de ]0, +∞[ dans ]0, 1[, d’inverse
1
G(u) = − ln(1 − u). (4.6)
λ
Lemme 4. Si U est de loi uniforme sur [0, 1], alors G(U ) suit la loi exponentielle de paramètre
λ.

En utilisant ce lemme on déduit l’algorithme suivant pour la simulation de la loi expo-

nentielle :

Algorithme Matlab :

function x=exponentielle(a)
% tirage suivant la loi exponentielle de parametre a
x=- log(rand())/a;
end

Ce procédé annonce la méthode de simulation par inversion de la fonction de répartition,

que nous présenterons plus tard.

Remarque 4. Pourquoi a-t-on remplacé 1 − U par U dans l’algorithme ?

4.2.6 Loi géométrique de paramètre p

Soit p ∈]0, 1]. On veut simuler une variable aléatoire géométrique de paramètre p, donc
de loi
+∞
X
µ(x) = (1 − p)k−1 pδk (x). (4.7)
k=1

On a donc dans ce cas les probabilités pk = P(X = k) = (1−p)k−1 p pour tout k ≥ 1 Plusieurs
méthodes sont possibles.

Méthode 1 : En utilisant la loi exponentielle :

h i
ln U
Lemme 5. Si U suit la loi uniforme sur [0, 1], alors X = 1 + ln(1−p) est une variable
aléatoire qui suit la loi géométrique de paramètre p.
h i
ln U
Démonstration : Notons X = 1 + ln(1−p) . Par définition de la partie entière, X prend ses

10
valeurs dans N∗ . Soit maintenant k ∈ N∗ . Calculons :
   
ln U
P(X = k) = P 1 + =k
ln(1 − p)
  
ln U
= P =k−1
ln(1 − p)
 
ln U
= P k−1≤ <k
ln(1 − p)
= P (k ln(1 − p) < ln U ≤ (k − 1) ln(1 − p))
= P (1 − p)k < U ≤ (1 − p)k−1

= (1 − p)k−1 − (1 − p)k
= (1 − p)k−1 p.


L’algorithme correspondant s’écrit :

Algorithme Matlab :
function x=geometrique2(p)
% tirage suivant la loi geometrique de parametre p
x=1+floor(log(rand())/log(1-p));
Méthode 2 : En utilisant la loi de Bernoulli :
Lemme 6. Si (Xi )i∈N∗ sont des variables aléatoires i.i.d. de loi de Bernoulli de paramètre
p, alors N = min{i : Xi = 1} est une variable aléatoire qui suit une loi géométrique de
paramètre p.
L’algorithme correspondant s’écrit.

function x=geometrique1(p)
% tirage suivant la loi geometrique de parametre p
x=1;
while rand()>p, x=x+1;
end

Méthode 3 : En utilisant la méthode générale pour la simulation des variables aléatoires

discrètes, voir 4.3.1.

4.2.7 Loi de Poisson de paramètre λ

Soit λ ∈ R+ . On veut simuler une variable aléatoire X ∼ P(λ), donc de loi donnée par
+∞
X λk
µ(x) = exp(−λ) δk (x). (4.8)
k=0
k!
k
On est ici dans la situation pk = P(X = k) = λk! e−λ , pour tout k ≥ 0.
Une variable aléatoire poissonienne ne prend pas ses valeurs dans un ensemble fini, mais
on peut utiliser la méthode générale pour les variables aléatoires discrètes présentée dans le

11
paragraphe 4.3.1. La méthode proposée ici annonce la méthode de la fonction inverse.

Méthode 1 : Cumul des durées exponentielles.

Une première méthode pour la simulation de variables aléatoires de Poisson est issue d’une
des propriétés des processus de Poisson, il s’agit du cumul de durées exponentielles.
On sait que si des événements surviennent à des dates séparées par des durées exponen-
tielles de paramètre λ, le nombre d’événements survenant en une unité de temps suit une loi
de Poisson de même paramètre.

Lemme 7. Soient (Xi )i∈N∗ des variables aléatoires i.i.d. exponentielles de paramètre λ. La
variable aléatoire définie par

0 si X1 ≥ 1
N= Pi
max{i ≥ 1 : j=1 Xj ≤ 1} sinon,

suit une loi de Poisson de paramètre λ.

Démonstration : Soit k ∈ N. On évalue :

P Pk+1 
k
P(N = k) = P j=1 Xj ≤ 1 < j=1 Xj
Z
= λk+1 exp(−λ(x1 + · · · + xk+1 ))1{x1 +···+xk ≤1}
Rk+1
+

Z ×1{xk+1 >1−(x1 +···+xk )} dx1 . . . dxk+1

= λk exp(−λ(x1 + · · · + xk ))1{x1 +···+xk ≤1}
Rk+
Z +∞ 
× λ exp(−λxk+1 )dxk+1 dx1 . . . dxk
Z 1−(x1 +···+xk )

= λ exp(−λ(x1 + · · · + xk ))1{x1 +···+xk ≤1} exp(−λ)

k
Rk+
× exp(λ(x 1 + · · · + xk ))dx1 . . . dxk
= λ exp(−λ) Rk 1{x1 +···+xk ≤1} dx1 . . . dxk .
k
R
+

Pour calculer la dernière intégrale, on effectue le changement de variable s1 = x1 , s2 =

x1 + x2 , . . . , sk = x1 + · · · + xk :
Z Z
1{x1 +···+xk ≤1} dx1 . . . dxk = 1{s1 ≤s2 ≤···≤sk ≤1} ds1 . . . dsk
Rk+ Rk+
Z 1 Z sk Z s2
1
= dsk dsk−1 . . . ds1 = .
0 0 0 k!
λk −λ
On voit donc que P(N = k) = k!
e . Ce qui termine la preuve. 

Remarque 5. Il faut donc simuler des variables aléatoires exponentielles de paramètre λ et

compter le nombre de simulations nécessaires pour dépasser 1, ou bien simuler des variables
aléatoires exponentielles de paramètre 1 et compter le nombre de simulations nécessaires pour
dépasser λ.

12
Remarque 6. Soient (Ui )i∈N∗ des variables aléatoires i.i.d. de loi uniforme sur [0, 1], alors
si Xi = − λ1 ln(Ui ), les (Xi )i∈N∗ sont des variables aléatoires i.i.d. exponentielles de paramètre
λ, et !
i i i
X 1 Y Y
Xj ≤ 1 ⇔ − ln Uj ≤ 1 ⇔ Uj ≥ exp(−λ).
j=1
λ j=1 j=1

L’algorithme s’écrit :

Algorithme Matlab :

function x=poisson(a)
% tirage suivant la loi de Poisson de parametre a
test=exp(-a);x=0;prod=rand();
while (prod>=test), x=x+1; prod=prod*rand(); end;

Méthode 2 : En utilisant la méthode générale pour les variables aléatoires discrètes, voir le
paragraphe 4.3.1.
On sait que

λk −λ
X ∼ P(λ) ⇔ pk = P{X = k} = e (pour k ∈ N)
k!
ce qui implique que
λ
pk+1 = pk ,
k+1
Pk
et en notant qk la somme des pj pour 0 ≤ j ≤ k, (qk = j=0 pj ), on a

λ
qk+1 = qk + pk .
k+1
On simule donc une variable de Poisson de paramètre λ en prenant

k 1{qk−1 <U ≤qk }

X
X=
k≥0

avec la convention q−1 = 0.

Cet algorithme est assez simple à programmer.

4.2.8 Loi gaussienne centrée réduite : Méthode de Box-Muller

La loi normale n’a pas une densité à support compact et on ne connaı̂t pas d’expression
simple de l’inverse de sa fonction de répartition. Nous utiliserons le résultat suivant pour sa
simulation :

Lemme 8. Soient U1 et U2 deux variables √ aléatoires de loi uniforme sur

√ [0, 1], supposées de
plus indépendantes. Alors, si on pose X = −2 ln U1 cos(2πU2 ) et Y = −2 ln U1 sin(2πU2 ),
les variables aléatoires X et Y sont i.i.d. de loi gaussienne centrée réduite.

Ceci conduit à une méthode de simulation simple pour une loi gaussienne, basée sur la
simulation de deux variables aléatoires uniformes sur [0, 1]. C’est la méthode de Box-Muller.

13
Algorithme Matlab :

function [x,y]=boxmuller
% tirage de deux N(0,1) independantes
r=sqrt(-2*log(rand())); t=2*pi*rand();
x=r*cos(t); y=r*sin(t);
end

4.3 Méthodes générales

4.3.1 Méthode générale pour une variable aléatoire discrète
Soit X une variable aléatoire sur un espace de probabilité (Ω, F, P), telle que X(Ω) =
{x0 , x1 , ...} = {xi }i∈I avec I = N ou I = {0, 1, . . . , n}. Pour tout i ∈ I, pi = P(X = xi ).
Considérons le cas I = {0, . . . , n} fini, le second se traitant de la même manière. La loi de
probabilité de X est donnée par
n
X
µ(x) = pi δxi (x)
i=0

où n
X
pi = P(X = xi ) = µ(xi ), ∀i ∈ I et pi = 1. (4.9)
i=0

Notons par qk le cumul des pi , pour 0 ≤ i ≤ k, i.e.

k
X
qk = pi . (4.10)
i=0

On a alors q0 = p0 et qn = 1.
Nous souhaitons simuler une variable aléatoire de même loi que X. Pour ce faire on va
construire une variable aléatoire Y telle que

P(X = xi ) = P(Y = xi ) = pi , ∀i ∈ I,

à l’aide d’une variable aléatoire U de loi uniforme sur [0, 1].

Algorithme L’algorithme s’écrit : on tire U une uniforme sur [0, 1] et on pose


Y = x0 si U ≤ p0 = q0
Y = xk si qk−1 < U ≤ qk , pour k ∈ {0, . . . , n − 1}

ce qui s’exprime également sous la forme

n
Y = x0 1{U ≤q0 } + xi 1{qi−1 <U ≤qi } .
X

i=1

Remarque 7. Cette méthode générale peut s’appliquer pour tous les variables aléatoires
discrètes particulières que nous avons étudiées.

14
 L’algorithme suivant simule une variable représentative pour cette loi :

Algorithme Matlab :
function y=simuldiscrete(x,p)
% simulation d’une v.a. de loi discrete
% x=vecteur des valeurs prises, p=vecteur des probabilites
u=rand(); q=p(1); i=0;
while (u>q); i=i+1; q=q+p(i); end;
y=x(i);
end
Démonstration : En effet, on a : P(X = x0 ) = P(rand() ≤ p0 ) = p0 et pour k ≥ 1,
k−1 k
!
X X
P(X = xk ) = P pi < rand() ≤ pi = pk .
i=0 i=0


Remarque 8. Le nombre N de tests nécessaires satisfait N = 1 ssi u ≤ p0 , et pour i > 1,
N = i ⇔ qi−1 < u ≤ qi .
On a donc intérêt à réordonner les (xi )i≥0 dans l’ordre des (pi )i≥0 décroissants.

4.3.2 Méthode de simulation par inversion de la fonction de répartition

On veut simuler une variable aléatoire X de fonction de répartition F .
La méthode utilisée pour simuler une loi exponentielle est en fait générale : dès que l’on
sait inverser une fonction de répartition F , il est très facile de simuler une variable aléatoire
de fonction de répartition F .
Lemme 9. Soit U une variable aléatoire qui suit la loi uniforme sur [0, 1], et F une fonction
de répartition bijective de ]a, b[ dans ]0, 1[ d’inverse F −1 . Alors F −1 (U ) est une variable
aléatoire de fonction de répartition F .
Démonstration : On pose X = F −1 (U ). Cette variable aléatoire prend ses valeurs dans
]a, b[. Remarquons que nécessairement F est strictement croissante de ]a, b[ dans ]0, 1[. Soit
t ∈]a, b[ :
P(X ≤ t) = P(F −1 (U ) ≤ t) = P(U ≤ F (t)) = F (t).
Donc la fonction de répartition de X est bien F . 
Remarque 9. L’hypothèse de la connaissance de F −1 n’a de sens que si F est strictement
croissante. Cependant, même dans ce cas, il se peut que F −1 n’ait pas d’expression analytique
simple, c’est le cas par exemple pour la loi normale.
Si on note φ(x) la fonction de répartition de la loi normale centrée réduite,
Z x
1 t2
φ(x) = √ e− 2 dt,
2π −∞
il n’existe pas de formulation simple de φ(x) et encore moins de φ−1 (x), la méthode de la
fonction inverse ne peut donc pas s’appliquer directement à la loi normale.
Il existe cependant des polynômes donnant de bonnes approximations de φ(x) et de φ−1 (x)
qui permettent donc d’appliquer la méthode de la fonction inverse à la loi normale moyennant
cette approximation.

15
4.3.3 Algorithme par rejet
On veut simuler une variable aléatoire X de densité f et de fonction de répartition F .
Exemple 3. Commençons par un exemple très simple : comment simuler une loi uniforme
sur le disque unité {x2 + y 2 ≤ 1} ?

function X=disque
% simule un point uniformement sur le disque unite
X=2*rand(1,2)-[1,1];
while (norm(X)>1),
X=2*rand(1,2)-[1,1];
end

L’idée est la suivante : on tire des points uniformément dans le carré [0, 1] × [0, 1], et on
les jette jusqu’à en obtenir un qui tombe dans le disque. La loi du point obtenue est la loi
d’un point tiré uniformément dans le carré conditionnellement à être dans le disque, ce qui
est encore la loi uniforme sur le disque.
Remarque 10. Quelle est la loi du nombre N de passages dans la boucle ?
Lemme 10. Soit X une variable aléatoire de densité f (sur R) à simuler. On suppose qu’il
existe une constante k > 0 et une densité g (sur R aussi, facile à simuler) tels que
∀x, f (x) ≤ kg(x).
Soit U une variable aléatoire de loi uniforme sur [0, 1] et Z une variable aléatoire, indépendante
de U , de densité g. On pose V = kU g(Z). Alors, la loi de Z conditionnellement à l’événement
{V < f (Z)} a pour densité f .
Remarque 11. Notons que nécessairement k ≥ 1 (car f, g sont des densités).
Remarque 12. Il est très important de noter qu’on doit choisir la constante k la plus petite
possible pour minimiser le nombre de rejets : plus la majoration est grossière, plus il faut des
tirages pour obtenir une valeur acceptable.
f (z)
Démonstration : Notons que pour tout z ∈ R, kg(z)
≤ 1. On a tout d’abord

P(V < f (Z)) = P(kU g(Z)

Z <1 f (Z)) 
Z
= g(z) 1{kug(z)<f (z)} du dz
ZR 0
f (z)
= g(z) dz
R kg(z)
1
= .
k
Évaluons ensuite
Z t Z 1 
P({Z ≤ t} ∩ {V < f (Z)}) = g(z) 1{kug(z)<f (z)} du dz
Z−∞
t
0
f (z)
= g(z) dz
−∞
Z t kg(z)
1
= f (z)dz.
k −∞

16
 Z t
Ce qui montre que P(Z ≤ t|V < f (Z)) = f (z)dz, donc la loi conditionnelle de Z sachant
−∞
que {V < f (Z)} a bien pour densité f .
Ce résutltat se généralise à Rd . 
On obtient donc l’algorithme de simulation par rejet (on suppose qu’on possède une fonc-
tion simulg qui simule une variable aléatoire de densité g) :

Algorithme :
function z=simulf
% simule par rejet une va de densite f
u=rand();
z=simulg;
v=k*u*g(z);
while (v>=f(z));
u=rand();
z=simulg;
v=k*u*g(z);
end
Démonstration : Notons N le nombre de tests fait lors de cette fonction. N est une va-
riable aléatoire, à valeurs dans N∗ . Notons (Un )n≥1 la suite des appels à la fonction rand(),
et (Zn )n≥1 la suite des appels à la fonction simulg. Toutes ces variables aléatoires sont
indépendantes, les premières de loi uniforme sur [0, 1], les secondes de densité g. On note
Vn = kUn g(Zn ), et on note X la sortie de la fonction.
Soit t ∈ R. Evaluons
1 t
Z
P(X ≤ t et N = 1) = P(V1 < f (Z1 ) et Z1 ≤ t) = f (z)dz
k −∞
par la démonstration précédente. Soit maintenant i ≥ 2. Par indépendance, et comme
précédemment :
P(X ≤ t et N = i)

= P(V1 ≥ f (Z1 ), V2 ≥ f (Z2 ), . . . , Vi−1 ≥ f (Zi−1 ), Vi < f (Zi ), Zi ≤ t)

= P(V1 ≥ f (Z1 ))P(V2 ≥ f (Z2 )) . . . P(Vi−1 ≥ f (Zi−1 ))P(Vi < f (Zi ), Zi ≤ t)

 i−1 Z t
1 1
= 1− f (z)dz.
k k −∞

Finalement (rappelons que k ≥ 1) :

+∞
X
P(X ≤ t) = P(X ≤ t et N = i)
i=1
+∞
X i−1 Z t
1 1
= 1− f (z)dz
i=1
k k −∞
Z t
= f (z)dz,
−∞
donc la densité de X est bien f . 

17
4.3.4 Simulation par composition
Exemple 4. Soit F et G deux fonctions de répartition sur R. On construit une variable
aléatoire X de la façon suivante : on lance une pièce qui tombe sur pile avec probabilité
1/3 et sur face avec probabilité 2/3, si pile sort, on tire un nombre au hasard suivant la loi
donnée par F , sinon, on tire un nombre au hasard suivant la loi donnée par G. Déterminer
la fonction de répartition H de X.

L’exemple précédent est un exemple de mélange de variables aléatoires. On suppose

Pn main-
tenant qu’on veut simuler une variable
P aléatoire X de fonction de répartition F = i=1 θi Fi ,
où les θi sont des poids : θi ≥ 0 et ni=1 θi = 1 et les Fi sont des fonctions de répartition dont
les lois sont faciles à simuler. On suppose qu’on a à notre disposition des fonctions simulFi
qui simulent des variables aléatoires de fonction de répartition Fi et une fonction simulTheta
qui simule une variable aléatoire Θ à valeur dans {1, ..., n} telle que P(Θ = i) = θi .

function x = melange
% simulation par melange
i=simulTheta;
x=simulFi;
end

Cette méthode se généralise immédiatement à un nombre infini dénombrable de poids

(θi )i∈N , et même à un mélange ”continu”
R de lois : on suppose qu’on veut simuler une variable
aléatoire X de densité f (z) = Θ g(θ)fθ (z)dθ, où g est une densité de probabilité, ainsi que
tous les fθ . L’algorithme est alors simple : on tire θ suivant g, puis x suivant fθ .

4.4 Simulation de vecteurs aléatoires

4.4.1 Cas indépendant
Supposons qu’on souhaite simuler Z un vecteur aléatoire de Rd . Si ses composantes
sont indépendantes, on est ramené au cas de la simulation de variables aléatoires réelles
indépendantes traité dans les sections précédentes (on rappelle que les sorties successives de
rand donnent des variables aléatoires indépendantes de loi uniforme sur [0, 1]). Le problème
est différent quand les coordonnées ne sont pas indépendantes.

4.4.2 Vecteur gaussien

Un vecteur gaussien dans Rd est caractérisé par un vecteur moyenne m ∈ Rd , et une
matrice de covariance K, de taille d × d, symétrique et positive. On suppose pour l’instant
que K est définie positive.

Théorème 1 (Décomposition de Cholesky). Si K est une matrice symétrique définie positive

de taille d × d, il existe au moins une matrice réelle triangulaire inférieure L telle que :

K = Lt L.

On peut également imposer que les éléments diagonaux de la matrice L soient tous strictement
positifs, et la factorisation correspondante est alors unique.

18
En pratique on cherche L par coefficients indéterminés et identification, colonne par colonne.
Voir l’exemple juste après.

Remarque 13. Si K est seulement semi-définie positive, la décomposition de Cholesky existe

encore mais elle n’est plus unique.

Lemme 11. Soit T un vecteur gaussien de dimension d, centré et de matrice de covariance

Id (dont les coordonnées sont des variables aléatoires i.i.d. normales centrées réduites). Soit
m ∈ Rd , et K une matrice définie positive de taille d × d. Soit L la matrice triangulaire
inférieure donnée par la factorisation de Cholesky de K.
Alors le vecteur aléatoire Z = m + LT est un vecteur gaussien de moyenne m et de
matrice de covariance K.

Démonstration : Comme L est une application linéaire de Rd dans Rd , Z est encore un

vecteur gaussien d-dimensionnel. Pour l’espérance :

E(Z) = E(m + LT ) = m + LE(T ) = m,

puisque T est centré. Pour la matrice de covariance, comme celle de T est Id ,

E((Z − m)t (Z − m)) = E(LT t T t L) = LE(T t T )t L = LIdt L = Lt L = K.

Exemple 5. On veut simuler le vecteur gaussien Z de R3 de moyenne m et de matrice de

covariance K avec    
1 1 −1 0
m =  −2  et K =  −1 5 6 .
4 0 6 10
On commence par chercher la matrice de factorisation de Cholesky L par coefficients
indéterminés et identification :
 
l11 0 0
L =  l21 l22 0  .
l31 l32 l33

On calcule LLt et on identifie, colonne par colonne :

2
1. l11 = k11 = 1 ⇒ l11 = 1.
2. l11 l21 = k21 = −1 ⇒ l21 = −1.
3. l11 l31 = k31 = 0 ⇒ l31 = 0.
2 2
4. l21 + l22 = k22 = 5 ⇒ l22 = 2.
5. l21 l31 + l22 l32 = k32 = 6 ⇒ l32 = 3.
2 2 2
6. l31 + l32 + l33 = k33 = 1 ⇒ l33 = 1.
Donc    
1 1 0 0
Z =  −2  +  −1 2 0  T,
4 0 3 1
où T est un vecteur gaussien de R3 dont les trois composantes sont i.i.d. centrées réduites.

19
4.4.3 Cas général
1. Le cas d’une v.a. discrète à valeurs dans Rd se traite comme celui d’une variable aléatoire
réelle discrète.
2. La méthode de rejet a été présentée dans la cas d’une variable aléatoire à valeurs dans Rd .
3. On peut encore utiliser les lois conditionnelles. Nous allons illustrer cette méthode, dite
méthode récurrente, par un exemple :

Exemple 6. On veut simuler un vecteur (X, Y ) de loi uniforme sur le triangle ABC avec
A = (0, 0), B = (1, 1) et C = (0, 1).

On commence par déterminer la densité de cette loi :

f (x, y) = 210≤x≤1 10≤y≤1 1x≤y .

On peut alors calculer la densité de X :

fX (x) = 2(1 − x)10≤x≤1 ,

puis la densité de la loi conditionnelle de Y sachant X :

1
fY |X (y|x) = 1x≤y≤1 .
1−x
Pour la simulation, on procède maintenant de la façon suivante : on simule X suivant sa
densité (en utilisant par exemple la méthode de la fonction de répartition), on obtient une
1
valeur x, puis on simule Y suivant la densité fY |X (y|x) = 1−x 1x≤y≤1 (on reconnaı̂t par exemple
une loi usuelle).

Remarque 14. Remarquons que la seconde étape ressemble beaucoup à la simulation par
mélange.

20
4.5 Exercices
Exercice 1. Soit X une variable aléatoire de loi de Bernoulli de paramètre p.
1. Proposer une fonction qui construit une réalisation de X.
2. Proposer une fonction qui renvoie un k-échantillon de cette variable aléatoire.
3. Tester ces résultats en Matlab pour différentes valeurs de p et k.
Exercice 2. Simuler une variable aléatoire X de loi binomiale de paramètres (n, p) pour
différentes valeurs de n et p.
Exercice 3. Soit X une variable aléatoire de loi uniforme sur {1, . . . n}. Proposer un algo-
rithme en Matlab pour la simulation de cette loi :
1. En utilisant la simulation d’une loi de probabilité discrète sur un ensemble fini.
2. En utilisant le résultat suivant : Si U est une variable aléatoire de loi uniforme sur [0, 1]
alors [nU ] + 1 est une variable aléatoire de loi uniforme sur {1, . . . , n}.
Simuler les résultats pour différentes valeurs de n. Comparer les deux approches.
Exercice 4. Soit X une variable aléatoire de loi uniforme sur [a, b] avec a, b ∈ R, a < b.
Proposer un algorithme pour la simulation de cette loi et tester le pour différentes valeurs de
a et b.
Exercice 5. Soit X une variable aléatoire de loi exponentielle de paramètre λ. Proposer un
algorithme pour la simulation de cette loi en utilisant la méthode de la fonction inverse.
Exercice 6. Soit X une variable aléatoire géométrique de paramètre p. Proposer un algo-
rithme pour la simulation de cette loi :
1. En utilisant le résultat : Si (Yi )i∈N∗ sont des variables aléatoires i.i.d. de loi de Bernoulli
de paramètre p, alors N = min{i; Yi = 1} suit une loi géométrique de paramètre p.
 
ln U
2. En utilisant le résultat : Si U suit une loi uniforme sur [0, 1] alors G = 1 +
ln(1 − p)
suit la loi géométrique de paramètre p.
3. Comparer les résultats en insistant sur le temps de calcul.
Exercice 7. Soit Y une variable aléatoire de loi exponentielle de paramètre λ, avec λ > 0.
Sa densité est donnée par
f (x) = λe−λx 1{x>0} .
On définit une variable aléatoire Z = 1 + [Y ], où, pour tout x ∈ R, [x] note la partie entière
de x.
1. Quelle est la loi de Z ? Justifier votre réponse.
2. Soient p ∈]0, 1[ et U une variable aléatoire de loi uniforme sur [0, 1]. Trouver la constante
λ telle que  
ln U
1+ − ∼ G(p).
λ
3. Utiliser ce résultat pour écrire un algorithme pour la simulation d’une variable aléatoire
X de loi géométrique de paramètre p.
4. Donner des résultats numériques de votre algorithme pour p = 0.1, p = 0.5 et p = 0.9
(4 valeurs pour chaque choix de p).

21
Exercice 8. Soit X une variable aléatoire de loi de Poisson de paramètre λ. Simuler cette
variable aléatoire en utilisant le résultat suivant :
Soit (Xi )i∈N∗ une suite de variables aléatoires i.i.d. de loi exponentielle de paramètre λ. Alors
la variable aléatoire :

0, si X1 ≥ 1
N= Pi
max{i ≥ 1; j=1 Xj ≤ 1}, sinon
suit une loi de Poisson de paramètre λ.
Exercice 9. Soit X une variable aléatoire de loi normale de moyenne m et variance σ 2 .
Simuler cette variable aléatoire en utilisant la méthode de Box-Muller.
Exercice 10. En utilisant la méthode de la fonction de répartition, simuler une loi de Cauchy.
Soit X une variable aléatoire de loi de Cauchy de paramètre a > 0. Sa densité est donée par :
a
f (x) = , x ∈ R.
π(x + a2 )
2

1. Vérfier que f est une densité.

2. Calculer sa fonction de répartition F .
(Indication : F (x) = π1 arctan( xa ) + 21 .)
3. En déduire G l’inverse de la fonction F .
4. Soit U une variable aléatoire de loi uniforme sur [0, 1]. En utilisant la périodicité de la
fonction tan, montrer que
  
1
tan π U − = tan(πU ), en loi .
2
(Indication : Montrer



que pour toute fonction ϕ continue et bornée on a
E ϕ tan π U − 12

= E{ϕ(tan(πU ))}.)
5. Écrire un algorithme pour simuler une variable aléatoire de loi de Cauchy de paramètre
a.
6. Donner des réalisations de votre code Matlab pour a = 10 et a = 1 (4 valeurs
numériques pour chaque choix de a).
7. Tracer la vraie densité et un histogramme à l’aide de cet algorithme pour a = 10 et
a = 1. Utiliser un échantillon de taille 10000.
Exercice 11. Simulation de la gaussienne par rejet par rapport à la double exponentielle.
On pose  2
1 z 1
f (z) = √ exp − et g(z) = exp(−|z|).
2π 2 2
1. Montrer que g est bien une densité sur R.
2. Déterminer une constante k satisfaisantq ∀z ∈ R, f (z) ≤ kg(z). On aura intérêt à
2e
prendre k la plus petite possible : [k = π
= 1, 3155].
3. Écrire un algorithme de simulation d’une loi gaussienne centrée réduite par rejet.
Exercice 12. Considérons l’algorithme de la méthode par rejet. Montrer que le nombre
de passage dans la boucle suit une loi géométrique dont le paramètre ne dépend que de k.
Justifier alors l’intérêt de prendre k la plus petite possible.

22
Exercice 13. Simuler un mélange d’exponentielles F (x) = α(1 − exp(−ax)) + (1 − α)(1 −
exp(−bx)), avec α ∈ [0, 1].

Exercice 14. Simuler une variable aléatoire de loi 21 δ0 (x) + 12 e−x 1R+ (x)dx.

Exercice 15. Simuler le vecteur gaussien Z de R3 de moyenne m et de matrice de covariance

K avec    
1 1 −1 0
m =  −2  et K =  −1 5 6 .
4 0 6 10

Exercice 16. Simuler un vecteur (X, Y ) de loi uniforme sur le triangle ABC avec A = (0, 0),
B = (1, 1) et C = (0, 1).

Exercice 17. Simuler un vecteur de loi f (x, y, z)dxdydz, où

f (x, y, z) = 61x>0,y>0,z>0 1x+y+z<1 :

1. En utilisant les lois conditionnelles.

2. Par la méthode du rejet.

3. Comparer.

Exercice 18. Simuler un vecteur de loi f (x, y)dxdy, où

1
f (x, y) = (x + y)e−(x+y) 1x>0 1y>0 .
2

23
 Simulation de variables aléatoires
Dans tout ce qui suit U note une v.a. de loi uniforme sur [0, 1] et (Un )n≥1 une suite de
v.a. i.i.d., avec U1 de loi uniforme sur [0, 1].

Loi Méthode de simulation

Bernoulli B(p), p ∈ [0, 1] 1{U ≤p}

n
1{Ui ≤p}
X
Binomiale B(n, p), n ∈ N∗ , p ∈ [0, 1]
i=1

Uniforme sur {0, . . . , n − 1} [nU ]

Uniforme sur [a, b] a + (b − a)U

ln(U )
Exponentielle E(λ) −
λ

 
ln(U )
Géométrique G(p) 1+
ln(1 − p)

( n
)
Y
Poisson P(λ) min n ∈ N : Ui ≤ e−λ
i=1

p
Gaussienne N (m, σ 2 ) m+σ −2 ln(U1 ) sin(2πU2 )

Cauchy C(a) a tan(πU )

24