Probabilités V

STEP, MINES ParisTech∗

10 janvier 2024 (#0258370)

Table des matières

Objectifs d’apprentissage 2

Intégrale de Monte-Carlo 3

Génération de nombres pseudo-aléatoires 4

Générateur de nombres uniformes pseudo-aléatoires . . . . . . . . 5

Méthodes de simulation de variables aléatoires réelles 6

Méthode d’inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Cas où FX est bijective . . . . . . . . . . . . . . . . . . . . . . . 6
Réciproque généralisée . . . . . . . . . . . . . . . . . . . . . . . . 8
Conséquences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Méthode d’inversion . . . . . . . . . . . . . . . . . . . . . . . . . 8
Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Méthode de rejet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
Loi uniforme sur un borélien . . . . . . . . . . . . . . . . . . . . . 9
Stabilité par conditionnement . . . . . . . . . . . . . . . . . . . . 9
Loi uniforme sur le sous-graphe d’une densité . . . . . . . . . . . 9
Méthode de rejet . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Simulation de variables aléatoires gaussiennes : Box-Muller . . . . . . 11
Proposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Simulation d’un vecteur gaussien à densité 12

∗ Ce document est un des produits du projet ‡ paulinebernard/CDIS issu de la collaboration

de (P)auline Bernard (CAS) et (T)homas Romary (GEOSCIENCES). Il dérive du projet ‡

boisgera/CDIS, initié par la collaboration de (S)ébastien Boisgérault (CAOR), (T)homas
Romary et (E)milie Chautru (GEOSCIENCES), (P)auline Bernard (CAS), avec la contribution
de Gabriel Stoltz (Ecole des Ponts ParisTech, CERMICS). Il est mis à disposition selon les
termes de la licence Creative Commons “attribution – pas d’utilisation commerciale – partage
dans les mêmes conditions” 4.0 internationale.

1
Echantillonnage d’importance 13
Echantillonnage d’importance . . . . . . . . . . . . . . . . . . . . 15
Densité optimale . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Exercices 16
Loi uniforme dans un domaine . . . . . . . . . . . . . . . . . . . . . . 16
Simulation selon la loi géométrique . . . . . . . . . . . . . . . . . . . . 16
Simulation de la loi gaussienne par la méthode de rejet . . . . . . . . . 18
Simulation selon la loi de Wigner . . . . . . . . . . . . . . . . . . . . . 18
Loi des grands nombres et théorème central limite . . . . . . . . . . . 18
Simulation d’un mélange de gaussiennes . . . . . . . . . . . . . . . . . 19
Echantillonnage d’importance . . . . . . . . . . . . . . . . . . . . . . . 19

Solutions 19
Loi uniforme dans un domaine . . . . . . . . . . . . . . . . . . . . . . 22
Simulation selon la loi géométrique . . . . . . . . . . . . . . . . . . . . 22
Simulation de la loi gaussienne par la méthode de rejet . . . . . . . . . 24
Simulation de la loi de Wigner . . . . . . . . . . . . . . . . . . . . . . 25
Loi des grands nombres et théorème central limite . . . . . . . . . . . 25
Simulation d’un mélange de gaussiennes . . . . . . . . . . . . . . . . . 25
Echantillonnage d’importance . . . . . . . . . . . . . . . . . . . . . . . 25

Annexe 25
Preuve de la méthode d’inversion . . . . . . . . . . . . . . . . . . . . . 25
Proposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Références 28

Objectifs d’apprentissage
Cette section s’efforce d’expliciter et de hiérarchiser les acquis d’apprentissages
associés au chapitre. Ces objectifs sont organisés en paliers :
(◦) Prérequis (•) Fondamental (••) Standard (•••) Avancé (••••) Expert
Sauf mention particulière, la connaissance des démonstrations du document n’est
pas exigible 1

Intégrale de Monte Carlo

— • connaître le principe de l’intégration par la méthode Monte Carlo
— • savoir que cette approche se justifie par la loi des grands nombres
— • savoir que l’approximation fournie par le TCL fournit un contrôle de
l’erreur
1. l’étude des démonstrations du cours peut toutefois contribuer à votre apprentissage, au
même titre que la résolution d’exercices.

2
Génération de nombres pseudo-aléatoires
— • connaître le principe de la génération de nombres pseudo-aléatoires par
la méthode des congruences

Méthodes de simulation de v.a.

— • connaître et savoir implémenter en python la méthode d’inversion
— • connaître et savoir implémenter en python la méthode de rejet
— • connaître la méthode de Box-Muller
— • connaître et savoir implémenter la simulation de vecteurs gaussien par
la méthode de Cholesky

Echantillonnage d’importance
— • connaître la définition de la méthode d’échantillonnage d’importance
— • savoir qu’un bon choix de densité instrumentale permet de minimiser
la variance d’estimation

Intégrale de Monte-Carlo
Les méthodes de Monte-Carlo ont été développées initialement par des physiciens
dans les années 1950 (notamment par les travaux de Metropolis, Ulam, Hastings,
Rosenbluth) pour calculer des intégrales (déterministes) à partir de méthodes
probabilistes numériquement assez économiques. Le nom a été donné en référence
au célèbre casino du fait du caractère aléatoire de ces méthodes.
Les méthodes de simulation sont basées sur la production de nombres aléatoires,
distribués selon une certaine loi de probabilité. Dans de nombreuses applications,
pour une certaine fonction h, on souhaite calculer, pour une variable aléatoire
X de loi PX Z
I = E (h(X)) = h(x) PX (dx).

En général, même si on sait évaluer h en tout point, on ne peut pas calculer for-
mellement l’intégrale I. Le calcul d’intégrale par la méthode Monte-Carlo consiste
dans sa version la plus simple à générer un échantillon (X1 , . . . , Xn ) ∼i.i.d. PX ,
et à estimer I par la moyenne empirique
n
1X
Mn (h) = h(Xi ),
n i=1

où i.i.d signifie indépendant et identiquement distribué. En effet, d’après la loi

forte des grands nombres, si h(x) est PX -intégrable, on a l’assurance que quand
n → +∞, Z
Mn (h) → h(x)PX (dx) p.s.

3
Si de plus h(X)2 est intégrable, la vitesse de convergence de Mn (h) peut être
évaluée, puisque la variance
Z
1 2
V(Mn (h)) = (h(x) − I) PX (dx)
n
peut également être estimée à partir de l’échantillon (X1 , . . . , Xn ) par la quantité
n
1X 2
σn2 = (h(Xi ) − Mn (h)) .
n i=1
Le théorème central limite assure alors que pour n grand,
Mn (h) − I
σn
suit approximativement une loi N (0, 1) 2 . Cette propriété conduit à la construc-
tion de tests de convergence et de bornes de confiance asymptotiques pour Mn (h).
Par exemple, on aura
P (I ∈ [Mn (h) − 1.96σn , Mn (h) + 1.96σn ]) ≈ 0.95,
−1
où 1.96 ≈ Φ (0, 975) avec Φ la fonction de répartition de la loi normale centrée
réduite.

√ propriété indique que la vitesse de convergence de Mn (h) est

En outre, cette
de l’ordre de n, et ce indépendamment de la dimension du problème. Cela
explique la supériorité de cette méthode par rapport aux méthodes d’intégration
numérique déterministes dont les vitesses de convergence décroissent rapidement
(exponentiellement) avec la dimension du problème.

Exercice – Calcul numérique d’une probabilité (••) Calculer une esti-

mation de la probabilité de l’événement sin(X) > 1/2 pour X de loi N (0, 1).
Indication : si x est un vecteur, y=([Link](x)>1/2) désigne le vecteur de
booléens dont les coordonnées valent True ou False selon que les coordonnées
correspondantes de [Link](x) soient ou non > 1/2. Par ailleurs, les vecteurs
de booléens se convertissent naturellement en vecteurs de 0 et de 1 si on leur
applique une fonction dont les arguments sont des nombres (ex : [Link]).
(Solution p. 19.)

Génération de nombres pseudo-aléatoires

Les ordinateurs sont des machines déterministes. Il peut sembler paradoxal de
leur demander de générer des nombres aléatoires 3 . En réalité, les algorithmes de
2. ce résultat sera démontré dans le cours de science des données au second semestre.
3. Von Neumann (1951) résume ce problème très clairement : “Any one who considers
arithmetical methods of reproducing random digits is, of course, in a state of sin. As has
been pointed out several times, there is no such thing as a random number — there are only
methods of producing random numbers, and a strict arithmetic procedure of course is not such
a method.”

4
génération de nombres aléatoires vont produire des séquences de nombres déter-
ministes qui vont avoir l’aspect de l’aléatoire. On s’intéresse ici à la génération
de nombres uniformes sur ]0, 1[ (on exclut les bornes par commodité), puisque,
comme on le verra dans la suite, il s’agit de l’ingrédient de base de toutes les
méthodes de simulation stochastique.

Définition – Générateur de nombres uniformes pseudo-aléatoires

Un générateur de nombres uniformes pseudo-aléatoires est un algorithme qui
étant donné une valeur initiale u0 et une transformation T produit une séquence
ui = T i (u0 ), i ∈ N∗ , de valeurs dans ]0, 1[.
Pour tout n ∈ N∗ , les valeurs (u1 , . . . , un ) reproduisent le comportement d’une
suite de variables aléatoires (V1 , . . . , Vn ) i.i.d de loi uniforme sur ]0,1[, lorsqu’on
les compare au travers d’un ensemble de tests statistiques 4 , vérifiant par exemple
que la corrélation entre deux nombres successifs est suffisamment faible.

Exemple – la méthode des congruences Cet algorithme, dû à Lehmer

(1951), est l’un des premiers à avoir été proposé et implémenté. Il repose sur 2
paramètres :
— le multiplicateur a ∈ N∗ ,
— le modulo m ∈ N∗ .
Etant donné un entier k0 tel que 1 ≤ k0 ≤ m − 1, la séquence de nombres est
générée par la transformation suivante : pour tout n ∈ N∗ on définit l’entier
kn+1 = a kn mod m puis on pose

kn+1 k0
un+1 = , u0 = .
m m

On peut remarquer que l’algorithme va produire une séquence de valeurs qui

sera périodique, c’est-à-dire qu’après un certain nombre d’itérations, la suite se
répétera, et qu’il pourra fournir au plus m − 1 valeurs différentes. Ce trait est
commun à tous les générateurs de nombre aléatoires et est lié aux limitations
matérielles des ordinateurs (on ne peut représenter qu’un nombre fini de nombres).
Le choix des valeurs de a et m est par conséquent crucial. Il existe des critères
qui permettent de s’assurer du bon comportement de cette suite :
— m est un nombre premier (le plus grand possible),
— p = (m − 1)/2 est un nombre premier,
— ap = −1 mod m.
Ce critère assure une période pleine et donc que tous les nombres (1/m, . . . , (m −
1)/m) seront générés. Cependant, les nombres générés ne sont pas indépendants,
pas même non corrélés. On peut montrer que la corrélation entre 2 nombres
successifs vaut approximativement 1/a. Il convient donc de choisir a suffisamment
4. Par exemple, la suite de tests Die Hard, due à Marsaglia.

5
grand pour que celle-ci devienne négligeable. Par exemple, en prenant a = 1000
et m = 2001179, on obtient une période de 2001178 et une corrélation de l’ordre
de 10−3 .
Ce générateur de nombres uniformes pseudo-aléatoires, bien que rudimentaire,
a ouvert la voie à des générateurs plus sophistiqués, toujours basés sur des
opérations arithmétiques. Le plus couramment utilisé à ce jour est l’algorithme
de Mersenne-Twister. Il s’agit du générateur par défaut de la plupart des logiciels
tels que Python, R, MATLAB, Julia, MS Excel,. . . Il passe notamment avec succès
toute la batterie de tests Die Hard. En particulier, sa période vaut 219937 − 1.
Pour certains usages, cet algorithme n’est cependant pas recommandé du fait de
sa prédictibilité. C’est notamment un défaut rédhibitoire pour les applications
en cryptographie. Il existe des variantes mieux adaptées à ce cas de figure. On
notera enfin qu’il existe des générateurs de nombres aléatoires basés sur des
phénomènes physiques comme un bruit électrique ou des phénomènes quantiques
et donc parfaitement imprévisibles.

Méthodes de simulation de variables aléatoires

réelles
On a vu au chapitre II du cours de Probabilités que l’on pouvait transformer
des variables aléatoires réelles suivant certaines lois pour en obtenir de nouvelles.
Par exemple, si X1 , . . . , Xn sont n ∈ N∗ variables gaussiennes centrées réduites
indépendantes, alors X12 + . . . , Xn2 suit une loi du χ2 à n degrés de liberté. Dans
le même esprit, on va voir ici comment simuler des v.a.r. de lois diverses à partir
de la simulation de variables uniformes sur ]0, 1[. On introduit une notation qui
sera utile dans la suite : pour spécifier que deux v.a.r. X et Y ont même loi, on
L
écrira X = Y .

Méthode d’inversion
L’objectif de ce paragraphe est de définir quand et comment il est possible de
simuler une variable aléatoire réelle X de fonction de répartition (f.d.r.) FX en
transformant la simulation d’une variable aléatoire U de loi uniforme sur ]0, 1[.
En d’autres termes, on cherche à déterminer les conditions sous lesquelles il est
L
possible d’identifier une fonction borélienne ψ : ]0, 1[→ R telle que X = ψ(U ).
Commençons par un cadre simple, où FX est bijective d’un intervalle non vide
de R sur ]0, 1[.

Proposition – Cas où FX est bijective

Soient X une variable aléatoire réelle de fonction de répartition FX et U une
variable uniforme sur ]0, 1[. S’il existe un intervalle non vide ]a, b[⊂ R tel que

6
 −1
FX : ]a, b[→ ]0, 1[ est bijective, de bijection réciproque FX : ]0, 1[→]a, b[, alors
−1 L L
FX (U ) = X et FX (X) = U .

Démonstration Le premier résultat est immédiat : pour tout x ∈ R, par

−1
croissance de FX et donc de FX , on a
−1


P FX (U ) ≤ x = P (U ≤ FX (x)) = FX (x).

Concernant le second, notons G la fonction de répartition de la variable aléatoire

FX (X). La croissance de FX garantit que P (FX (X) ∈ ]0, 1[) = 1. Ainsi, pour
tout x ∈ R on a bien
1 si x ≥ 1,
−1


G(x) = P (FX (X) ≤ x) = P X ≤ FX (x) = x si 0 < x < 1,
0 sinon.

Exercice – Exemples d’application 1 (•) Donner un algorithme de simula-

tion d’une v.a.r. X suivant une loi
— uniforme sur un intervalle I ⊂ R,
— Exponentielle de paramètre λ ∈ R∗+ ,


−1
— de Cauchy, de densité x ∈ R 7→ π 1 + x2 ,
∗
— de Laplace
n de paramètres
o µ ∈ R et s ∈ R + , de densité x ∈ R 7→
1
exp − |x−µ|
2s s ,
— Logistique de paramètres µ ∈ R et s ∈ R∗+ , de fonction de répartition

−1
x ∈ R 7→ 1 + exp − x−µ

s ?
(Solution p. 20.)
−1
Dans cette situation idéale, ψ = FX est une solution à notre problème. Que se
passe-t-il en revanche si FX n’est pas bijective ?

Exercice – Pile ou face (•) Proposer une méthode pour simuler un tir à
pile ou face à partir de la simulation d’une variable uniforme sur ]0, 1[. (Solution
p. 20.)
On a déjà vu au chapitre II que les fonctions de répartition de v.a.r. possèdent
un nombre au plus dénombrable de points de discontinuité. Sur chaque intervalle
où elles sont continues, on peut alors considérer qu’elles sont bijectives, quitte à
réduire les zones de palier à un point. Cela permet de généraliser la notion de
bijection réciproque pour ces fonctions.

7
Définition – Réciproque généralisée
Soit F une fonction de répartition. On définit sa réciproque généralisée (aussi
appelée inverse généralisée ou pseudo-inverse) comme la fonction

F − : u ∈ ]0, 1[ 7→ inf {x ∈ R : F (x) ≥ u} ∈ R.

Remarque – Conséquences
— Cette fonction est bien définie sur tout ]0, 1[, car quel que soit u dans cet
intervalle, l’ensemble {x ∈ R : F (x) ≥ u} est non vide et minoré. S’il était
vide ou non minoré pour un certain u0 ∈ ]0, 1[, pour tout x ∈ R on aurait
dans le premier cas F (x) < u0 < 0 et dans le second F (x) ≥ u0 > 1.
L’une comme l’autre de ces inégalités est impossible pour une fonction de
répartition.
— La réciproque généralisée de la f.d.r. FX d’une v.a.r. X est aussi
appelée
− 1
fonction quantile. On pourra notamment remarquer que FX 2 est la
médiane de X.
— Lorsque F réalise une bijection d’un intervalle non vide I ⊂ R sur ]0, 1[,
sa réciproque généralisée coïncide avec sa bijection réciproque.
−
On a alors le résultat suivant, qui stipule que ψ = FX est une solution universelle
à notre problème. La preuve détaillée est donnée en Annexe.

Théorème – Méthode d’inversion

Soient U une variable uniforme sur ]0, 1[ ainsi que X une variable aléatoire
−
réelle de fonction de répartition FX et de réciproque généralisée FX . Alors
− L
FX (U ) = X.

Exercice – Exemples d’application 2 (•) Donner un algorithme de simula-

tion d’une v.a.r. X suivant une loi
— Binomiale de paramètres n ∈ N∗ et p ∈ ]0, 1[,
— Uniforme sur l’union de deux segments non vides et disjoints [a, b], [c, d] ⊂
R, de densité x ∈ R 7→ (b − a + d − c)−1 1[a,b]∪[c,d] (x) ?
(Solution p. 20.)

Limitations
La méthode d’inversion peut sembler universelle pour simuler une v.a.r. X à
partir de U ∼ U]0,1[ . Cependant, elle nécessite en pratique de disposer d’une
expression analytique de FX pour en déduire sa réciproque généralisée. Ce n’est
typiquement pas le cas de nombreuses lois usuelles comme la loi Normale. On va
donc déterminer d’autres procédures pour simuler des variables suivant de telles
lois.

8
Méthode de rejet
La méthode de rejet est une alternative populaire à la méthode d’inversion,
lorsque cette dernière ne peut être utilisée directement et que la loi cible
possède une densité. On la doit à Von Neumann (1951). Pour en comprendre
le fondement, il nous faut d’abord introduire une généralisation naturelle de la
loi uniforme dans R à tout Rd (d ∈ N∗ ). On notera λ la mesure de Lebesgue sur
Rd .

Définition – Loi uniforme sur un borélien

La loi uniforme sur un borélien A ⊂ Rd de volume λ(A) > 0 est une loi de
probabilité admettant pour densité

1A (x)
f : x ∈ Rd 7→ .
λ(A)

Exercice – Loi uniforme sur un pavé (•) Comment simuler un vecteur

aléatoire (U1 , . . . , Ud ) de loi uniforme sur un pavé non vide [a1 , b1 ]×· · ·×[ad , bd ] ⊂
Rd ? (Solution p. 21.)
Une propriété fondamentale de cette loi est qu’elle reste stable par conditionne-
ment, dans le sens suivant.

Proposition – Stabilité par conditionnement

Soit U un vecteur aléatoire de loi uniforme sur un borélien A ⊂ Rd de vo-
lume λ(A) > 0. Alors pour tout borélien B ⊂ A de volume λ(B) > 0, la loi
conditionnelle PU |U ∈B de U sachant que U ∈ B est uniforme sur B.

Exercice – Démonstration (•) (Solution p. 21.)

Exercice – Simulation (•) Déduire de la proposition précédente (p. 9) une

méthode pour simuler un vecteur aléatoire de loi uniforme sur un borélien B ⊂ Rd
borné et de volume λ(B) > 0. (Solution p. 21.)
Il est aussi possible de simuler un vecteur aléatoire de loi uniforme sur certains
boréliens non vides et non bornés, comme l’établit la proposition ci-dessous.

Proposition – Loi uniforme sur le sous-graphe d’une densité

Soient une densité f : R → R, une v.a.r. X et une variable U uniforme sur ]0, 1[,
indépendante de X. On note Af := {(x, y) ∈ R × R+ : f (x) ≥ y} le domaine
limité par le graphe de f et l’axe des abscisses (souvent appelé sous-graphe de f ).
Si X est de densité f , alors le couple (X, U f (X)) suit une loi uniforme sur Af .

9
Démonstration Commençons par remarquer que λ(Af ) = 1. Quel que soit
(z, v) ∈ R2 , par indépendance de X et U et par Fubini on a
Z Z
P (X ≤ z, U f (X) ≤ v) = 1]−∞,z] (x) 1]−∞,v] (uf (x)) 1]0,1[ (u) f (x) du dx
ZRz RZ 
= 1]−∞,v]∩]0,f (x)[ (uf (x)) f (x) du dx
−∞ R
Z z Z v
= 1]0,f (x)[ (u) du dx
−∞ −∞
Z z Z v
= 1Af (x, u) du dx.
−∞ −∞

Ainsi, (X, U f (X)) admet pour densité 1Af , qui correspond bien à celle d’une loi
uniforme sur Af . ■
Pour simuler un vecteur uniforme (X, Y ) sur un ensemble Af tel que défini à la
proposition précédente, il suffit donc de simuler une v.a.r. X de densité f , puis
une variable U uniforme sur ]0, 1[, et de poser Y = f (X)U .
En combinant les résultats précédents, on obtient la méthode de rejet, illustrée
sur la figure ci-dessous, où AY := {(x, y) ∈ R × R+ : y ≤ fY (x)}.

Méthode de rejet
On souhaite simuler une variable aléatoire réelle X de densité fX . Supposons
que l’on sait simuler une variable U uniforme sur ]0, 1[, ainsi qu’une v.a.r. Y
(par exemple avec la méthode d’inversion) de densité fY telle qu’il existe un
réel a > 0 pour lequel on a ∀x ∈ R : fX (x) ≤ a fY (x). Il suffit alors de suivre
l’algorithme suivant :

1. simuler Y et U ,
2. si aU fY (Y ) > fX (Y ), recommencer à l’étape 1.
3. poser X = Y .

Limitations
La méthode de rejet a l’avantage de permettre de simuler des variables aléatoires
à densité dont la fonction de répartition n’a pas de forme analytique, rendant la
méthode d’inversion inapplicable. Néanmoins, pour pouvoir l’appliquer il faut
connaître une densité auxiliaire qui, multipliée par un réel positif, majore la
densité cible, et que l’on sait simuler. Le taux de rejet, c’est-à-dire la probabilité
de l’événement {aU fY (Y ) > fX (Y )}, peut parfois être élevé, notamment lorsque
la dimension de X est grande, ce qui limite l’efficacité de la méthode.

Exercice – Taux de rejet (••) Calculer le taux de rejet de la méthode

proposée ci-dessus. (Solution p. 22.)

10
 y
AY
a fY (x)

(Y, aU fY (Y ))
fX (x)

x
0

Figure 1 – Méthode de rejet

Simulation de variables aléatoires gaussiennes : Box-Muller

On a vu que la méthode d’inversion est inappropriée pour simuler une variable
gaussienne, puisqu’elle requiert l’expression analytique de la fonction de répar-
tition cible. Il existe des méthodes basées sur une intégration numérique de la
densité gaussienne puis une inversion de cette approximation de la f.d.r. mais
elle ne sont pas optimales en temps de calcul. La méthode de rejet est quant
à elle sous-optimale, dans le sens où toutes les variables uniformes générées ne
sont pas directement utilisées (une partie, potentiellement grande, est rejetée).
La loi normale étant fondamentale en probabilité, il est plus que souhaitable de
pouvoir en trouver une méthode de simulation exacte et efficace.
George E. P. Box et Mervin E. Muller ont proposé en 1958 une telle méthode
(G.E.P. Box and M.E. Muller (1958)). Elle exploite la propriété d’invariance
par rotation de la densité d’un couple de variables gaussiennes indépendantes
centrées réduites.

Proposition – Proposition
Soient U et V deux variables
p indépendantes, de loi uniforme
p sur ]0, 1[. Alors les
variables aléatoires X = −2 ln(U ) cos (2πV ) et Y = −2 ln(U ) sin (2πV ) sont
indépendantes et suivent toutes deux une loi normale centrée réduite.

Démonstration On considère (X, e Ye ) un vecteur aléatoire dont les deux

composantes sont indépendantes, de loi normale centrée réduite. On note ses
coordonnées polaires aléatoires R e et Θ.
e On a vu dans le cours de Probabi-
lités II que dans ce cas, R e et Θe sont indépendantes, la première de densité
2
−r
fRe : r ∈ R 7→ r e 2 1R∗+ (r) et la seconde de loi uniforme sur ]0, 2π].
X et Y sont les coordonnées
p cartésiennes obtenues à partir d’un rayon et d’un
angle : en posant R = −2 ln(U ) et Θ = 2πV on obtient X = R cos(Θ) et
Y = R sin(Θ). Par indépendance de U et V , R et Θ sont indépendantes.
 Pour que

le vecteur (X, Y ) ait la même distribution que (X,
e Ye ) = R e R
e cos(Θ), e sin(Θ)
e ,

11
 L e L e
il suffit de montrer que R = R et Θ = Θ.
— Commençons par étudier la loi de R, depfonction de répartition FR .
On remarque que la fonction u ∈ ]0, 1[7→ −2 ln(u) ∈ R∗+ est bijective,
strictement décroissante. Ainsi, pour tout r ∈ R− on a P (R ≤ r) = 0 et
pour tout r ∈ R∗+ on a
r2 r2
p   
P(R ≤ r) = P −2 ln(U ) ≤ r = P U ≥ e− 2 = 1 − e− 2 .

qui correspond à la fonction de répartition associée à la densité fR

e.
Z r ir
x2 x2 r2
h
x e− 2 dx = −e− 2 = 1 − e− 2 si r > 0,
Z r
fR
e(x) dx = 0 0
−∞
0 sinon.

— En ce qui concerne la loi de Θ, de fonction de répartition FΘ , puisque la

fonction v ∈ ]0, 1[ 7→ 2πv ∈ ]0, 2π[ est bijective strictement croissante, on
a que pour tout θ ∈ R

1 si θ ≥ 2π,
θ

θ
FΘ (θ) = P V ≤ 2π = si θ ∈ ]0, 2π[,
2π
0 si θ ≤ 0,

qui est bien la fonction de répartition d’une loi uniforme sur ]0, 2π[.
■
Cette méthode permet de simuler directement deux variables gaussiennes centrées
réduites indépendantes à partir de deux variables uniformes indépendantes. Pour
simuler une variable gaussienne d’espérance m ∈ R et de variance σ 2 ∈ R∗+
quelconques, on se rappelera que si X suit une loi normale centrée réduite, alors
σX + m suit une loi normale d’espérance m et de variance σ 2 .

Simulation d’un vecteur gaussien à densité

On souhaite simuler un vecteur gaussien X = (X1 , . . . , Xd ) à valeurs dans Rd
d’espérance m et de matrice de covariance C définie positive (et donc inversible)
données.
Puisque la matrice C est définie positive, elle admet une racine carrée, c’est-à-dire
qu’il existe une matrice N telle que C = N N t . En effet, on peut par exemple
décomposer C de la manière suivante :

C = V D V t,

où V est une matrice orthogonale et D est la matrice diagonale dont les termes
diagonaux sont les valeurs propres (toutes strictement positives) de C. Il suffit

12
alors de prendre N = V D1/2 , où D1/2 est la matrice diagonale dont les termes
diagonaux sont les racines carrées des valeurs propres.
En pratique, il est coûteux numériquement d’effectuer le calcul des valeurs
propres et des vecteurs propres de C. On va plutôt calculer sa décomposition ou
factorisation de Cholesky qui permet d’écrire

C = L Lt

avec L une matrice triangulaire inférieure. 5

Soit maintenant un autre vecteur gaussien Y = (Y1 , . . . , Yd ) à valeurs dans Rd
et de matrice de covariance identité, notée Id . Autrement dit, les Yi sont des
variables aléatoires gaussiennes centrées, réduites et indépendantes.
Alors, le vecteur Z = m + L Y est gaussien, d’espérance m et de matrice de
covariance C. En effet, Z est gaussien comme combinaison linéaire de variables
aléatoires gaussiennes, E(Z) = E(m+L Y ) = m et V(Z) = E((L Y )2 ) = LId Lt =
C.

Echantillonnage d’importance
On introduit dans cette section la méthode d’échantillonnage d’importance
(importance sampling en anglais), que l’on appelle aussi, de manière plus intuitive,
échantillonnage préférentiel, pour les lois à densité. Pour ce faire, nous allons
commencer par un exemple qui montre qu’il peut être plus efficace de simuler
des valeurs selon une loi différente de celle d’intérêt, autrement dit de modifier
la représentation de l’intégrale I sous la forme d’une espérance calculée selon
une autre densité.

Exemple – Loi de Cauchy conditionnelle Supposons que l’on s’intéresse à

calculer la probabilité p qu’une variable X de loi de Cauchy standard soit plus
grande que 2 (on peut le calculer directement et p = 0.15)
Z +∞
1
p= dx.
2 π(1 + x2 )

Si on estime p directement à partir d’un échantillon (X1 , . . . , Xn ) simulé selon

la loi de Cauchy standard, soit
n
1X
pb1 = 1X >2 ,
n i=1 i

5. Cette décomposition est très utile dans la résolution de systèmes linéaires de la forme
A x = b, où b est connu, x inconnu et A est définie positive. Cela revient à résoudre L Lt x = b.
On pose alors y = Lt x et on résout d’abord Ly = b, ce qui est très rapide puisque L est
triangulaire inférieure (on commence par y1 = b1 /L11 , puis y2 = (b2 − L21 y1 )/L22 , etc. en
descendant). On résout ensuite Lt x = y, ce qui est aussi très rapide pour la même raison (on
commence par xn = yn /Lnn puis on remonte).

13
la variance de l’estimateur V(bp1 ) = p(1 − p)/n = 0.127/n, puisque pb1 suit une loi
binomiale de paramètre (n, p). On peut réduire cette variance (et donc améliorer
la qualité de l’estimateur) en tirant parti de la symétrie de la densité de la loi de
Cauchy, en formant un second estimateur
n
1X
pb2 = 1|Xi |>2 ,
n i=1

p2 ) = p(1 − p/2)/2n = 0.052/n. La relative inefficacité

dont la variance vaut V(b
de ces méthodes est due au fait que la majeure partie des valeurs simulées seront
en dehors de la zone d’intérêt ]2, +∞[. En passant par le complémentaire, on
peut réécrire p comme
Z 2
1 1
p= − dx,
2 0 π(1 + x2 )
2
dont le second terme peut être vu comme l’espérance de h(U ) = π(1+U 2 ) avec

U ∼ U[0,2] . Tirant un échantillon (U1 , . . . , Un ) i.i.d. de loi uniforme sur [0, 2], on
obtient un troisième estimateur :
n
1 1X 2
pb3 = − ,
2 n i=1 π(1 + Ui2 )

dont la variance vaut V(b p3 ) = (E(h(X)2 ) − E(h(U ))2 )/n = 0.0285/n (par
intégration par parties). Enfin, on peut encore réécrire (voir Ripley (1987))
1/2
y −2
Z
p= dy,
0 π(1 + y −2 )
 
V −2
qui peut être vue comme E 2π(1+V −2 ) avec V ∼ U]0,1/2[ . L’estimateur formé
à partir de cette représentation et d’un échantillon (V1 , . . . , Vn ) i.i.d. de loi
de 0.95 10−4 /n. Il est donc bien plus efficace
uniforme sur [0, 1/2] a une variance √
que pb1 puisqu’il nécessite environ 103 ≈ 32 fois moins de simulations pour
atteindre la même précision.
On a ainsi vu sur ce cas particulier que l’estimation d’une intégrale de la forme
Z
I = E (h(X)) = h(x)f (x)dx,
Rd

peut s’écrire de différentes manières, en faisant varier h et f . Par conséquent, un

estimateur “optimal” devrait tenir compte de l’ensemble de ces possibilités. C’est
justement l’idée développée dans la méthode d’échantillonnage d’importance
dont le principe est décrit dans la définition suivante.

14
Définition – Echantillonnage d’importance
La méthode d’échantillonnage d’importance est une évaluation de I basée sur la
simulation d’un échantillon X1 , . . . , Xn de loi de densité g et approximant :
n
1 X f (Xi )
Ef (h(X)) ≈ h(Xi ),
n i=1 g(Xi )

où la notation Ef signifie que l’espérance est calculée avec X ∼ f . On appelle

souvent les ratios fg(X
(Xi )
i)
les poids d’importance que l’on note wi .
Cette méthode est basée sur la représentation suivante de I :
Z
f (x)
Ef [h(X)] = h(x) g(x) dx
g(x)

que l’on appelle l’identité fondamentale de l’échantillonnage d’importance et

l’estimateur converge du fait de la loi forte des grands nombres.
Cette identité indique qu’une intégrale du type I n’est pas intrinsèquement
associée à une loi donnée. L’intérêt de l’échantillonnage d’importance repose sur
le fait qu’il n’y a aucune restriction sur le choix de la densité g, dite instrumentale,
que l’on peut donc choisir parmi les densités des lois que l’on sait simuler aisément.
Il y a bien évidemment des choix qui sont meilleurs que d’autres. Remarquons
tout d’abord que bien que l’estimateur proposé dans la définition ci-dessus (p.
15) converge presque sûrement, sa variance est finie si

f 2 (X) f 2 (x)
    Z
f (X)
Eg h2 (X) 2 = Ef h2 (X) = h2 (x) dx < ∞.
g (X) g(X) g(x)

On préconise alors l’usage de densités instrumentales g dont la queue de distri-

bution est plus épaisse que celle de f pour éviter que cette variance puisse être
infinie (on notera que cela dépend aussi de la fonction h à intégrer). En pratique,
on utilise généralement l’estimateur suivant, de variance finie et qui donne des
résultats plus stables numériquement que celui de la définition :
Pn
w h(Xi )
Pn i
i=1
i=1 wi

où on a remplacé n par la somme des poids d’importance. Puisque n1 i=1n wi =

P
1
Pn f (xi )
n i=1 g(xi ) tend vers 1 quand n → ∞, cet estimateur converge presque sûre-
ment vers Ef (h(X)) par la loi forte des grands nombres (voir C. P. Robert and
G. Casella (2004)).
Parmi les densités g qui fournissent des estimateurs de variance finie, il est
possible d’exhiber la densité optimale (au sens de la variance de l’estimateur
associé) pour une fonction h et une densité f données.

15
Théorème – Densité optimale
Le choix de g qui minimise la variance de l’estimateur donné dans la définition
ci-dessus (p. 15) est
|h(x)|f (x)
g ∗ (x) = R .
|h(x)|f (x)dx

Démonstration Notons d’abord que

2
f 2 (X)
    
h(X)f (X) 2 f (X)
Vg = Eg h (X) 2 − Eg h(X)
g(X) g (X) g(X)
et le second terme ne dépend pas de g. On minimise donc le premier terme.
D’après l’inégalité de Jensen, on a
2 Z 2
f 2 (X)
  
f (X)
Eg h2 (X) 2 ≥ Eg |h(X)| = |h(x)|f (x)dx ,
g (X) g(X)
qui nous donne une borne inférieure indépendante du choix de g. Elle est atteinte
en prenant g = g ∗ . ■
Ce résultatR formel n’a que peu d’intérêt pratique : le choix optimal de g fait
intervenir |h(x)|f (x)dx, qui est à une valeur absolue près la quantité que l’on
souhaite estimer ! Il suggère néanmoins de considérer des densités g telles que
|h|f /g est quasi constante et de variance finie. On se reportera au chapitre 3 de C.
P. Robert and G. Casella (2004) pour des exemples où un bon choix de g permet
des améliorations considérables par rapport à des estimateurs de Monte-Carlo
plus naïfs.

Exercices
Loi uniforme dans un domaine
Ecrire un algorithme pour simuler un point uniforme dans les trois domaines
représentés ci-dessous.

Simulation selon la loi géométrique

Cette loi correspond à la loi du nombre d’essais avant d’obtenir un évènement
de probabilité p
X ∼ G(p) P(X = n) = p(1 − p)n−1 n = 1, 2, . . .

Question 1 Calculer l’espérance et la variance de X (Solution p. 22.)

Question 2 Générer un échantillon de taille 1000 pour p = 0.6, en utilisant

une approche génétique. Calculer la moyenne et la variance. Comparer avec les
valeurs théoriques. (Solution p. 24.)

16
 (0,1)

(0,0) (1,0)
Figure 2 – Domaine A

0.00

0.75
0.25

1.00
1.25
1.00
0.75

0.00

Figure 3 – Domaine B

(5,2)

(3,1)

(0,0) (4,0)

Figure 4 – Domaine C

17
Question 3 Montrer que si X ∼ E(λ), alors ⌈X⌉ ∼ G(p). On précisera la
valeur de λ. (Solution p. 24.)

Question 4 Donner, implémenter et tester un algorithme pour simuler G(p)

directement (Solution p. 24.)

Simulation de la loi gaussienne par la méthode de rejet

X est une variable aléatoire de loi gaussienne d’espérance m = 0 et de variance
σ 2 = 1, N (0, 1). Sa densité est
 
1 1
f (x) = √ exp − x2 , x ∈ R
2π 2

e
p
Question 1 Montrer que f (x) ≤ C exp (−|x|) où C = 2π (Solution p. 24.)

Question 2 On considère g(x) = 12 exp (−|x|) , x ∈ R. Montrer que g est une

densité de probabilité. (Solution p. 25.)

Question 3 Proposer un algorithme pour simuler selon g. (Solution p. 25.)

Question 4 Implémenter l’algorithme de rejet pour simuler selon f . (Solution

p. 25.)

Simulation selon la loi de Wigner

La loi de Wigner
√ (ou du demi-cercle) est la loi de support [−2, 2] et de densité
1
f (x) = 2π 4 − x2 . Simuler un grand nombre de variables aléatoires de loi de
Wigner, représenter l’histogramme associé et le comparer à la densité. (Solution
p. 25.)

Loi des grands nombres et théorème central limite

Question 1 Afin d’illustrer la Loi des Grands Nombres, visualiser la suite
Sn = X1 +...+X
n
n
pour (Xi )i∈N∗ une suite de variables aléatoires indépendantes
de loi uniforme sur [-1, 1].
Indication : pour x = [x1 , . . . , xn ] vecteur, [Link](x) est le vecteur [x1 , x1 +
x2 , x1 + x2 + x3 , . . . , x1 + . . . + xn ] des sommes cumulées des coordonnées de x.
(Solution p. 25.)

Question 2 Faire de même avec des variables aléatoires Yi de loi de Cauchy

1 Y1 +...+Yn
(c’est-à-dire de densité π(1+x 2 ) ). La suite n semble-t-elle converger ?
Pourquoi ? (Solution p. 25.)

18
Question 3 Calculer la moyenne µ et l’écart-type
√ σ des Xi et vérifier, afin
d’illustrer le Théorème Central Limite, que n(Sn −µ)/σ converge en loi, lorsque
n → ∞, vers une variable de loi N (0,
√ 1) : on se fixera une grande valeur de n, on
simulera un grand nombre de fois n(Sn − µ)/σ, on en tracera l’histogramme
2
et on le comparera à la densité (2π)− 1/2e−x la loi N (0, 1). (Solution p. 25.)

Simulation d’un mélange de gaussiennes

On souhaite simuler selon la loi de mélange introduite dans l’exercice Mélange
de loi du chapitre 3.

Question 1 Soit K une v.a.r. de loi :


 1/2
 si x = 1
1/4 si x = 2

P(K = x) =

 1/4 si x = 3
0 sinon


Implémenter la méthode d’inversion pour simuler K (Solution p. 25.)

Question 2 Soit X1 , X2 et X3 des variables aléatoires de lois respectives

N (0, 1), N (5, 1/2) et N (8, 4). Implémenter un algorithme de simulation de XK .
Comparer l’histogramme obtenu pour un échantillon de taille 1000 avec la densité
de XK . (Solution p. 25.)

Echantillonnage d’importance
On cherche à évaluer l’espérance d’une variable aléatoire X gaussienne centrée
réduite (de densité fX et de f.d.r. FX ) sachant qu’elle dépasse la valeur 3.

Question 1 Exprimer la densité conditionnelle de X|X > 3. (Solution p. 25.)

Question 2 Calculer la quantité souhaitée à partir d’un échantillon de taille

1000 généré selon FX . Quelle est la valeur du taux de rejet ? (Solution p. 25.)

Question 3 Implémenter l’algorithme d’échantillonnage d’importance en pre-

nant comme loi instrumentale la loi gaussienne d’espérance 3 et de variance 1.
Quelle est la proportion de poids nuls ?

(Solution p. 25.)

Solutions
Calcul numérique d’une probabilité cf. notebook

19
Exemples d’application 1 On suppose que U ∼ U]0,1[
— Uniforme sur un intervalle I ⊂ R :
Si I =]a, b[, a < b, alors la fonction de répartition de X ∼ UI est
FX (x) = x−a
b−a 1]a,b[ (x) + 1[b,+∞[ (x), d’où par inversion X = a + (b − a)U
— Exponentielle de paramètre λ ∈ R∗+ :
la fonction de répartition de X ∼ E(λ) est FX (x) = 1−exp(−λx), d’où par
inversion X = − ln(1−U λ
)
. Pour simplifier, on remarquera que si U ∼ U]0,1[
, alors 1 − U ∼ U]0,1[ , d’où X = − ln(U λ
)


−1
— de Cauchy, de densité x ∈ R 7→ π 1 + x2 ,
la loi de Cauchy admet la fonction de répartition F (x) = π1 arctan(x) + 12 ,
d’où X = tan(πU − 1/2)
— de Laplace
n de oparamètres µ ∈ R et s ∈ R∗+ , de densité x ∈ R 7→
1
2s exp − |x−µ|
s ,
la loi de Laplace admet la fonction de répartition
 1
2 exp(x) si x < 0, 1
F (x) = = (1 + sgn(x)(1 − exp(−|x|)))
1 − 12 exp(−x) si x ≥ 0, 2
d’où X = sgn(U − 1/2) ln(1 − 2|U − 1/2|), où sgn est la fonction signe.
— Logistique de paramètres µ ∈ R et s ∈ R∗+ , de fonction de répartition

−1
x ∈ R 7→ 1 + exp − x−µ

s :

par inversion, X = s − ln U1 − 1 + µ

Pile ou face Il suffit d’associer à “pile” un événement portant sur U ∼ U]0,1[

de la bonne probabilité, p ∈]0, 1[, par exemple {U ≤ p}. On obtient ainsi
l’algorithme :

1. Générer u selon U]0,1[

2. Renvoyer “pile” si u ≤ p, “face” sinon

Exemples d’application 2
— Binomiale de paramètres n ∈ N∗ et p ∈ ]0, 1[,
La fonction de répartition de la loi binomiale est donnée par

0P
 si x < 0
⌊x⌋ n
k n−k
F (x) = k=0 k p (1 − p) si 0 ≤ x < n

1 si x ≥ n


où ⌊x⌋ est la partie entière de x. Pm

On peut alors découper ]0, 1[ en intervalles de la forme ] k=0 nk pk (1 −


Pm+1

p)n−k , k=0 nk pk (1 − p)n−k [ pour m ∈ {0, n − 1} et renvoyer m + 1 si
un u généré selon U]0,1[ appartient à l’un de ces intervalles, 0 sinon.
Il est cependant beaucoup plus simple de remarquer qu’une variable
aléatoire de loi B(n, p) s’obtient comme la somme de n variables de
Bernoulli que l’on peut simuler comme ci-dessus (p. 20).

20
 — Uniforme sur l’union de deux segments non vides et disjoints [a, b], [c, d] ⊂
R, de densité x ∈ R 7→ (b − a + d − c)−1 1[a,b]∪[c,d] (x)
La fonction de répartition d’une telle variable est
x−a b−a x−c+b−a
F (x) = 1]a,b[ (x)+ 1[b,c] (x)+ 1]c,d[ (x)+1[d,+∞[ (x)
b−a+d−c b−a+d−c b−a+d−c
Pour simuler cette variable, on peut donc utiliser l’algorithme suivant :
1. Générer u selon U]0,1[
2. Retourner
(
b−a
a + u(b − a + d − c) si u ∈]0, b−a+d−c [
c + (u(b − a + d − c) − (b − a)) sinon

Loi uniforme sur un pavé On note [a1 , b1 ] × · · · × [ad , bd ] ⊂ Rd . On

peut voir que la densité d’une v.a. uniforme sur un tel pavé s’écrit f (x) =
1 1
b1 −a1 · · · bd −ad 1[a1 ,b1 ]×···×[ad ,bd ] (x) soit comme le produit des densités de ses
coordonnées qui sont donc indépendantes. On peut donc simuler la variable
d’intérêt avec cet algorithme :

1. Générer u1 , . . . , ud indépendamment selon U]0,1[

2. Retourner (a1 + u1 (b1 − a1 ), . . . , ad + ud (bd − ad ))

Démonstration Soit C ∈ B(B) (considérer les boréliens de B suffit à carac-

tériser la loi conditionnelle à l’événement U ∈ B puisque si C ∩ B = ∅, alors
nécessairement, P(U ∈ C|U ∈ B) = 0). On a

P(U ∈ C, U ∈ B) P(U ∈ C)
P(U ∈ C|U ∈ B) = =
P(U ∈ B) P(U ∈ B)
R dx
C λ(A)
=R dx
puisque U ∼ UA
B λ(A)
Z
λ(C) 1
= = dx
λ(B) C λ(B)

on reconnaît bien la loi uniforme sur B.

Simulation On prend un pavé R tel que B ⊂ R puis on applique l’algorithme

suivant :

1. générer u selon UR
2. Si u ∈
/ R, retour en 1.
3. retourner u

21
Taux de rejet Puisque le couple (Y, U fY (Y )) est uniforme sur Af , alors le
couple (Y, aU fY (Y )) est uniforme sur Aaf , d’où

P({aU fY (Y ) > fX (Y )}) = E(1{aU fY (Y )>fX (Y )} )

= E(E(1{aU fY (y)>fX (y)} |Y = y))
Z Z afY (y)
1
= 1{aufY (Y )>fX (Y )} dufY (y)dy
R 0 afY (y)
Z  
fX (y)
= 1− fY (y)dy
R afY (y)
1
=1−
a

Loi uniforme dans un domaine

cf. notebook

Simulation selon la loi géométrique

Question 1 P(X = n) = p(1−p)n−1 = (1−q)q n−1 , où p est la probabilité de succès, n =
1, 2, . . .
∞
X
E[X] = iP(X = i)
i=1
X∞
= i(1 − q)q i−1
i=1
∞
X
= (1 − q) iq i−1
i=1
∞
X d i
= (1 − q) q
i=1
dq
∞
!
d X
i
= (1 − q) q
dq i=1
  ∞ ∞
d q X X q
= (1 − q) , car qi = q qi =
dq 1−q i=1 i=0
1−q
1
= (1 − q)
(1 − q)2
1
=
(1 − q)
1
=
p

Pour la variance V[X] = E[X 2 ] − E[X]

22
 ∞
X
E[X 2 ] = i2 P(X = i)
i=1
∞
X
= i2 (1 − q)q i−1
i=1
∞
X
= (1 − q) i2 q i−1
i=1
∞
X
i (i + 1)q i−1 − q i−1



= (1 − q)
i=1
∞
X
i(i + 1)q i−1 − iq i−1


= (1 − q)
i=1
∞  2 
X d i+1 d i
= (1 − q) q − q
i=1
dq 2 dq
∞ ∞
! !!
d2 X i+1 d X
i
= (1 − q) 2
q − q
dq i=1
dq i=1

On a déjà calculé le second terme ci-dessus.

∞ ∞
d2 X i+1 d2 X i
q = q q
dq 2 i=1 dq 2 i=1
d2 q 2
=
dq 2 1 − q
Calcul des dérivées :
d q2 2q(1 − q) + q 2
=
dq 1 − q (1 − q)2
2q − q 2
=
(1 − q)2

d 2q − q 2 2(1 − q)(1 − q)2 + 2(1 − q)(2q − q 2 )

2
=
dq (1 − q) (1 − q)4
2(1 − q)2 + 2(2q − q 2 )
=
(1 − q)3
2 − 4q + 2q 2 + 4q − 2q 2
=
(1 − q)3
2
=
(1 − q)3

23
On en déduit
 
2 1
E[X 2 ] = (1 − q) 3
−
(1 − q) (1 − q)2
2 1
= −
(1 − q)2 (1 − q)
2 − (1 − q)
=
(1 − q)2
1+q
=
(1 − q)2
Et finalement

V[X] = E[X 2 ] − E[X]

1+q 1
= −
(1 − q)2 (1 − q)2
q
=
(1 − q)2
1−p
=
p2

Question 2 cf. notebook

Question 3

P(⌈X⌉ = n) = P(n − 1 < X ≤ n)

= FX (n) − FX (n − 1)
= e−λ(n−1) (1 − e−λ )

d’où λ = − ln(1 − p)

Question 4 cf. notebook

Simulation de la loi gaussienne par la méthode de rejet

X est une variable aléatoire de loi gaussienne d’espérance m = 0 et de variance
σ 2 = 1, N (0, 1). Sa densité est
 
1 1 2
f (x) = √ exp − x , x ∈ R
2π 2

Question 1 Soit x ∈ R

x2 ≥ 2|x| − 1
r
1 −x2 /2 e −|x|
√ e ≤ e
2π 2π

24
Question 2 g est bien positive et son intégrale sur R vaut
Z Z Z
1 −|x|
g(x)dx = e dx = e−x dx = 1
R R 2 R+

Question 3 voir dans exemples d’application 1 (p. 7).

Question 4 cf. notebook

Simulation de la loi de Wigner

cf. notebook

Loi des grands nombres et théorème central limite

Question 1 cf. notebook

Question 2 cf. notebook

Question 3 cf. notebook

Simulation d’un mélange de gaussiennes

Question 1 cf. notebook

Question 2 cf. notebook

Echantillonnage d’importance
P(X≤x,X>3)
Question 1 On a FX|X>3 (x) = P(X ≤ x|X > 3) = 1−FX (3) =
FX (x)−FX (3) fX (x)
1−FX (3) 1]3,+∞[ (x) d’où fX|X>3 (x) = 1−FX (3) 1]3,+∞[ (x)

Question 2 cf. notebook

Question 3 cf. notebook

Annexe
Preuve de la méthode d’inversion
Pour pouvoir démontrer le théorème de la méthode d’inversion (p. 8), il faut
d’abord établir un certain nombre de propriétés de la réciproque généralisée
d’une fonction de répartition. Elle peuvent être visualisées sur la figure ci-dessous.

25
Proposition – Proposition
Soit F une fonction de répartition. Alors sa réciproque généralisée F − satisfait
les propriétés suivantes.

1. F − est croissante.
2. ∀ x ∈ R : F − ◦ F (x) ≤ x.
3. ∀ u ∈ ]0, 1[ : F ◦ F − (u) ≥ u avec égalité si u ∈ F (R).
4. ∀ (u, x) ∈ ]0, 1[ ×R : {F (x) ≥ u} ⇔ {x ≥ F − (u)} et {F (x) < u} ⇒
{x ≤ F − (u)}.

Démonstration Pour tout u ∈ ]0, 1[ on note F −1 ([u, 1]) := {x ∈ R : F (x) ≥ u}

l’image réciproque de [u, 1] par F .

1. Soit (u, v) ∈ ]0, 1[2 . Si u < v alors F −1 ([v, 1]) ⊂ F −1 ([u, 1]) d’où F − (u) =
inf F −1 ([u, 1]) ≤ inf F −1 ([v, 1]) = F − (v). La fonction F − est donc bien
croissante.
2. Soit x ∈ R, alors F − ◦ F (x) = inf z ∈ R : F (z) ≥ F (x) =


inf F −1 ([F (x), 1]). Comme F est croissante sur R on a [x, +∞[ ⊆
F −1 ([F (x), 1]) donc F − ◦ F (x) = inf F −1 ([F (x), 1]) ≤ inf[x, +∞[ = x.
3. Soit u ∈ ]0, 1[.
— Puisque F −1 ([u, 1]) est nécessairement non vide, il existe une suite décrois-
sante (xn )n∈N ⊆ F −1 ([u, 1]) convergeant vers inf
F −1 ([u, 1]) = F − (u).
La croissance de F implique que la suite F (xn ) n∈N est elle aussi dé-
croissante, minorée par u car (xn )n∈N ⊆ F −1 ([u, 1]) donc convergente.
Sa limite est de même supérieure ou égale à u. Comme F est continue à
droite, cette dernière n’est autre que
 
lim F (xn ) = F lim xn = F ◦ F − (u).
n→+∞ n→+∞

Nous avons donc bien F ◦ F − (u) ≥ u.

— Supposons maintenant que u ∈ F (R). Alors F −1 ({u}) := {x ∈ R : F (x) = u} =
̸
∅. Il existe donc une suite décroissante (xn )n∈N ⊆ F −1 ({u}) convergeant
vers inf F −1 ({u}) = F − (u) par croissance de F . Comme F est continue
à droite en tout point de R et F (xn ) = u pour tout n ∈ N, on a bien
 
u = lim F (xn ) = F lim xn = F ◦ F − (u).
n→+∞ n→+∞

4. Soit (u, x) ∈ ]0, 1[ ×R.

— Equivalence. Supposons F (x) ≥ u. On a F − ◦ F (x) ≥ F − (u) par
croissance de F − (propriété 1.) et F − ◦ F (x) ≤ x d’après la propriété 2.
Réciproquement, supposons que x ≥ F − (u). Alors par croissance de F
on a F (x) ≥ F ◦ F − (u) puis F ◦ F − (u) ≥ u d’après la propriété 3.

26
 — Implication. Supposons F (x) < u. Tout z ∈ F −1 ([u, 1]) vérifie F (z) ≥
u > F (x). Comme F est croissante sur R on a donc z ≥ x, ce qui implique
x ≤ inf F −1 ([u, 1]) = F − (y).
■

x
=
y
2

)
(x
F−
1.75

1.5

1.25

1
F (x)

0.75

0.5

0.25

x
0.25 0.5 0.75 1 1.25 1.5 1.75 2

Figure 5 – Réciproque généralisée d’une fonction de répartition

Nous pouvons maintenant établir la preuve du théorème de la méthode d’inversion

(p. 8).

Démonstration – Méthode d’inversion Soient x ∈ R et (Ω, A, P) l’espace

probabilisé sur lequel sont définies U et X. D’après la propriété 4 ci-dessus, on a
−

ω ∈ Ω : FX (U (ω)) ≤ x = {ω ∈ Ω : U (ω) ≤ F (x)}, d’où
−


P FX (U ) ≤ x = P (U ≤ FX (x)) = FX (x).

27
Références
C. P. Robert, and G. Casella. 2004. Monte-Carlo Statistical Methods. 2nd edition,
Berlin: Springer.
G.E.P. Box, and M.E. Muller. 1958. “A Note on the Generation of Random
Normal Deviates.” Ann. Math. Stat. 29: 610–11.
Lehmer, Derrick H. 1951. “Mathematical Methods in Large-Scale Computing
Units.” Annu. Comput. Lab. Harvard Univ. 26: 141–46.
Ripley, Brian D. 1987. Stochastic Simulation. New York: Wiley.
Von Neumann, John. 1951. “Various Techniques Used in Connection with Ran-
dom Digits.” U.S. Nat. Bur. Stand. Math. Ser. 12: 36–38.

28