1

IFT-3655, Modèles Stochastiques

aft
Simulation: Méthode de Monte Carlo

Prof. Pierre L’Ecuyer

DIRO, Université de Montréal

Dr
Ces “diapos” sont un support pour les présentations en classe.
Un traitement beaucoup plus détaillé de ce sujet se trouve ici:
[Link]
 2

Simulation d’un chaı̂ne de Markov

Supposons que l’on veut simuler sur ordinateur l’évolution d’une chaı̂ne de Markov

aft
{Xn , n ≥ 0} avec espace d’états X = {1, . . . , k} et probabilités de transition Pi,j .
À chaque étape n, on est dans un état Xn−1 = i, et on veut générer le prochain état Xn selon
les probabilités Pi,j = P(Xn = j | Xn−1 = i). Comment faire cela?

Dr
 2

Simulation d’un chaı̂ne de Markov

Supposons que l’on veut simuler sur ordinateur l’évolution d’une chaı̂ne de Markov

aft
{Xn , n ≥ 0} avec espace d’états X = {1, . . . , k} et probabilités de transition Pi,j .
À chaque étape n, on est dans un état Xn−1 = i, et on veut générer le prochain état Xn selon
les probabilités Pi,j = P(Xn = j | Xn−1 = i). Comment faire cela?

Supposons par exemple que Pi,1 = 1/4, Pi,2 = 1/2, Pi,3 = 0, et Pi,4 = Pi,5 = 1/8.
Les logiciels pour la simulation possèdent des fonctions pour générer U ∼ U(0, 1).
Pour générer Xn , on peut générer U ∼ U(0, 1), et retourner Xn = j pour le j qui correspond
à l’intervalle où se trouve U:

0
Xn = 1
Dr 1
4
Xn = 2
3
4
Ces probabilités sont en général différentes selon l’état courant i.
Xn = 4
7
8
Xn = 5
1
 3

Problème du collectionneur: simulation

Il y a k types d’items et on tire un item à la fois, chacun étant du type j avec probabilité pj ,

aft
pour j = 1, . . . , k.
On peut vouloir simuler ce système n fois pour estimer par exemple la distribution du nombre
N de tirages requis pour obtenir une collection complète.
Pour simuler les tirages, on doit simuler une suite de réalisations indépendantes d’une variable
aléatoire Y telle que P(Y = j) = pj . Ce Y represente le prochain type d’item obtenu.
Pour tirer Y , on divise l’intervalle (0, 1) en sous-intervalles de longueurs p1 , p2 , . . . , pk , on
tire U ∼ U(0, 1), et on retourne le numéro j de l’intervalle correspondant.

0 p1
Dr
Posons Fj = P(Y ≤ j) = p1 + · · · + pj .

Y =1
F1
Y =2
p2
F2
Y =3
p3
F3
etc.
1
 4

Un réseau d’activités stochastique.

Le graphe donne des relations de précédence entre les activités.

aft
L’activité j a une durée aléatoire Yj (la longueur de l’arc j) de fonction de répartition (cdf)
Fj , i.e., P[Yj ≤ y ] = Fj (y ). Supposons pour simplifier que les Yj sont indépendants.
La durée T du projet est la longueur (aléatoire) du plus long chemin de la source au puit.

5 8 puit
Y10
Y5
Y9 Y12

Dr
source 0
Y1

Y0
Y2
2

1
Y3
Y4
4

3
Y7
Y8

Y6
7

6
Y11
 5

On peut vouloir estimer P[T > x] pour x fixé, ou encore toute la loi de probabilité de T .
Calcul exact? Trop difficile en général.

aft
Dr
 5

On peut vouloir estimer P[T > x] pour x fixé, ou encore toute la loi de probabilité de T .
Calcul exact? Trop difficile en général.

aft
Simulation: répéter n fois:
générer les Yj selon Fj pour chaque j, puis calculer T .
On va voir plus loin comment générer les Yj .
On obtient ainsi n réalisations indépendantes T1 , . . . , Tn de T .
On peut estimer la densité de T par un histogramme des Ti ,
puis E[T ] et P[T > x] par les moyennes empiriques, etc.

Dr
 5

On peut vouloir estimer P[T > x] pour x fixé, ou encore toute la loi de probabilité de T .
Calcul exact? Trop difficile en général.

aft
Simulation: répéter n fois:
générer les Yj selon Fj pour chaque j, puis calculer T .
On va voir plus loin comment générer les Yj .
On obtient ainsi n réalisations indépendantes T1 , . . . , Tn de T .
On peut estimer la densité de T par un histogramme des Ti ,
puis E[T ] et P[T > x] par les moyennes empiriques, etc.

Illustration numerique:
Dr
Yj ∼ max(0, N(µj , σj2 )) pour j = 0, 1, 3, 10, 11, et Yj ∼ Expon(1/µj ) sinon.
µ0 , . . . , µ12 : 13.0, 5.5, 7.0, 5.2, 16.5, 14.7, 10.3, 6.0, 4.0, 20.0, 3.2, 3.2, 16.5.
On peut vouloir estimer P[T > 90] parce qu’on paye une pénalité si ça arrive.
 6
Idée naive: remplacer chaque Yj par son espérance. Donne T = 48.2.

aft
Dr
 6
Idée naive: remplacer chaque Yj par son espérance. Donne T = 48.2.
J’ai simulé de modèle n = 100 000 fois, en utilisant la librairie Java SSJ.
L’histogramme des n réalisations de T donne pas mal d’information.

aft
Les valeurs de T vont de 14.4 à 268.6; 11.57% dépassent x = 90.

fréquence (×103 )
mean = 64.2
T = 48.2
T = x = 90
10

0
0
Dr
25 50 75 100
ξˆ0.99 = 131.8

125 150 175 200

T
 7

Simulation de variables aléatoires uniformes indépendantes

Objectif: On veut produire des suites de valeurs qui ont l’air d’être tirées au hasard.

aft
Dr
 7

Simulation de variables aléatoires uniformes indépendantes

Objectif: On veut produire des suites de valeurs qui ont l’air d’être tirées au hasard.

aft
Exemple: Suites de bits (pile ou face):

011110100110110101001101100101000111?...
Loi uniforme: chaque bit est 1 avec probabilité 1/2.

Dr
 7

Simulation de variables aléatoires uniformes indépendantes

Objectif: On veut produire des suites de valeurs qui ont l’air d’être tirées au hasard.

aft
Exemple: Suites de bits (pile ou face):

01111?100110?1?101001101100101000111...
Loi uniforme: chaque bit est 1 avec probabilité 1/2.
Uniformité et indépendance:
Exemple: on 8 possibilités pour les 3 bits ? ? ?:

Dr 000, 001, 010, 011, 100, 101, 110, 111

On veut une proba. de 1/8 pour chacune, peu importe les autres bits.
 7

Simulation de variables aléatoires uniformes indépendantes

Objectif: On veut produire des suites de valeurs qui ont l’air d’être tirées au hasard.

aft
Exemple: Suites de bits (pile ou face):

01111?100110?1?101001101100101000111...
Loi uniforme: chaque bit est 1 avec probabilité 1/2.
Uniformité et indépendance:
Exemple: on 8 possibilités pour les 3 bits ? ? ?:

Dr 000, 001, 010, 011, 100, 101, 110, 111

On veut une proba. de 1/8 pour chacune, peu importe les autres bits.

Pour s bits, on veut une probabilité de 1/2s pour chacune des 2s possibilités.
 8

Loi uniforme sur (0, 1)

Pour la simulation en général, on veut imiter une suite U0 , U1 , U2 , . . . de variables aléatoires

aft
indépendantes de loi uniforme sur (0, 1).
Uniformité: On veut P[a ≤ Uj ≤ b] = b − a.

0 a b 1

Dr
 8

Loi uniforme sur (0, 1)

Pour la simulation en général, on veut imiter une suite U0 , U1 , U2 , . . . de variables aléatoires

aft
indépendantes de loi uniforme sur (0, 1).
Uniformité: On veut P[a ≤ Uj ≤ b] = b − a.

0 a b 1
Indépendance: Pour (U1 , . . . , Us ) en s dimensions, on veut
P[aj ≤ Uj ≤ bj pour j = 1, . . . , s] = (b1 − a1 ) · · · (bs − as ).

Exemple pour s = 2: 1

Dr U2
b2
a2

0 a1 b1 1 U 1
 8

Loi uniforme sur (0, 1)

Pour la simulation en général, on veut imiter une suite U0 , U1 , U2 , . . . de variables aléatoires

aft
indépendantes de loi uniforme sur (0, 1).
Uniformité: On veut P[a ≤ Uj ≤ b] = b − a.

0 a b 1
Indépendance: Pour (U1 , . . . , Us ) en s dimensions, on veut
P[aj ≤ Uj ≤ bj pour j = 1, . . . , s] = (b1 − a1 ) · · · (bs − as ).

Exemple pour s = 2: 1

Dr U2
b2
a2

0 a1 b1 1 U 1
Cette notion de v.a. uniformes et indépendantes n’est qu’une abstraction mathématique.
N’existe peut-être même pas dans la réalité physique!
 9

Générer un entier au hasard dans {0, 1, . . . , k − 1}

Générer U ∼ U(0, 1) et retourner X = ⌊kU⌋.

aft
Exemple avec k = 5:
0 1 2 3 4
0 1 2 3 4 1
k k k k

0 1 2 3 4
0

Dr
1 2

Valeur de V = kU.
3 4 5
 9

Générer un entier au hasard dans {0, 1, . . . , k − 1}

Générer U ∼ U(0, 1) et retourner X = ⌊kU⌋.

aft
Exemple avec k = 5:
0 1 2 3 4
0 1 2 3 4 1
k k k k

0 1 2 3 4
0

Dr
1 2

Valeur de V = kU.
3

Si on veut X uniforme sur {1, . . . , k}, on retourne X = 1 + ⌊kU⌋.

4 5
 10

Permutation
: aléatoire
1234567

aft
Dr
 10

Permutation
: aléatoire
1234567

aft
123467 5

Dr
 10

Permutation
: aléatoire
1234567

aft
123467 5
13467 52

Dr
 10

Permutation
: aléatoire
1 2 3 4567

aft
1 2 3 467 5
1 3 4 67 52
3 4 6 7 521

Dr
 10

Permutation
: aléatoire
1 2 3 4567

aft
1 2 3 467 5
1 3 4 67 52
3 4 6 7 521

Pour n objets, on choisit un entier de 1 à n,

puis un autre entier de 1 à n − 1, puis de 1 à n − 2, ...
On veut que chaque permutation ait la même probabilité.
Ex.: pour permuter 52 cartes, il y a 52! ≈ 2226 possibilités.

Dr
 11

Bits aléatoires par des mécanismes physiques

aft
Trajectoires de photons (système vendu par id-Quantique):

Dr
 12

Bruit thermique dans les résistances de circuits électroniques

aft
Dr temps
 12

Bruit thermique dans les résistances de circuits électroniques

aft
Dr0 1 0 1 0 0 1 1 1 0 0 1

On échantillonne le signal (voltage) périodiquement.

temps
 12

Bruit thermique dans les résistances de circuits électroniques

aft
Dr
00010110010100110 ···

On échantillonne le signal (voltage) périodiquement.

temps
 13

Plusieurs mécanismes sont brevetés et disponibles commercialement.

Aucun n’est parfait.

aft
Dr
 13

Plusieurs mécanismes sont brevetés et disponibles commercialement.

Aucun n’est parfait. On peut dimimuer le biais et/ou la dépendance en combinant des blocs

aft
de bits. Par exemple par un XOR:
0 1 |{z}
|{z} 1 0 |{z}
0 0 |{z}
1 0 |{z}
0 1 |{z}
1 0 |{z}
1 1 |{z}
0 1 |{z}
00
1 1 0 1 1 1 0 1 0

Dr
 13

Plusieurs mécanismes sont brevetés et disponibles commercialement.

Aucun n’est parfait. On peut dimimuer le biais et/ou la dépendance en combinant des blocs

aft
de bits. Par exemple par un XOR:
0 1 |{z}
|{z} 1 0 |{z}
0 0 |{z}
1 0 |{z}
0 1 |{z}
1 0 |{z}
1 1 |{z}
0 1 |{z}
00
1 1 0 1 1 1 0 1 0

ou encore (élimine le biais):

0 1 |{z}
|{z} 1 0 |{z}
0 0 |{z}
1 0 |{z}
0 1 |{z}
1 0 |{z}
1 1 |{z}
0 1 |{z}
00
0 1 1 0 1 0

Dr
 13

Plusieurs mécanismes sont brevetés et disponibles commercialement.

Aucun n’est parfait. On peut dimimuer le biais et/ou la dépendance en combinant des blocs

aft
de bits. Par exemple par un XOR:
0 1 |{z}
|{z} 1 0 |{z}
0 0 |{z}
1 0 |{z}
0 1 |{z}
1 0 |{z}
1 1 |{z}
0 1 |{z}
00
1 1 0 1 1 1 0 1 0

ou encore (élimine le biais):

0 1 |{z}
|{z} 1 0 |{z}
0 0 |{z}
1 0 |{z}
0 1 |{z}
1 0 |{z}
1 1 |{z}
0 1 |{z}
00
0 1 1 0 1 0

Mécanisme physique:
Dr
Essentiel pour cryptologie, loteries, etc. Mais pas pour la simulation.
Encombrant, pas reproduisible, pas toujours fiable, pas d’analyse mathématique de
l’uniformité et de l’indépendance à long terme.
 14

Générateurs algorithmiques (pseudo-aléatoires, GPA)

Mini-exemple: On veut imiter des nombres de 1 à 100 tirés au hasard.

aft
Dr
 14

Générateurs algorithmiques (pseudo-aléatoires, GPA)

Mini-exemple: On veut imiter des nombres de 1 à 100 tirés au hasard.

aft
1. Choisir un nombre x0 au hasard dans {1, . . . , 100}.
2. Pour n = 1, 2, 3, ..., retourner xn = 12 xn−1 mod 101 .

Dr
 14

Générateurs algorithmiques (pseudo-aléatoires, GPA)

Mini-exemple: On veut imiter des nombres de 1 à 100 tirés au hasard.

aft
1. Choisir un nombre x0 au hasard dans {1, . . . , 100}.
2. Pour n = 1, 2, 3, ..., retourner xn = 12 xn−1 mod 101 . Par exemple, si x0 = 1:

x1 = (12 × 1 mod 101) = 12,

Dr
 14

Générateurs algorithmiques (pseudo-aléatoires, GPA)

Mini-exemple: On veut imiter des nombres de 1 à 100 tirés au hasard.

aft
1. Choisir un nombre x0 au hasard dans {1, . . . , 100}.
2. Pour n = 1, 2, 3, ..., retourner xn = 12 xn−1 mod 101 . Par exemple, si x0 = 1:

x1 = (12 × 1 mod 101) = 12,

x2 = (12 × 12 mod 101) = (144 mod 101) = 43,

Dr
 14

Générateurs algorithmiques (pseudo-aléatoires, GPA)

Mini-exemple: On veut imiter des nombres de 1 à 100 tirés au hasard.

aft
1. Choisir un nombre x0 au hasard dans {1, . . . , 100}.
2. Pour n = 1, 2, 3, ..., retourner xn = 12 xn−1 mod 101 . Par exemple, si x0 = 1:

x1 = (12 × 1 mod 101) = 12,

x2 = (12 × 12 mod 101) = (144 mod 101) = 43,
x3 = (12 × 43 mod 101) = (516 mod 101) = 11, etc.
xn = 12n mod 101.

Visite tous les nombres de 1 à 100 une fois chacun avant de revenir à x0 .

Dr
(Parce que m = 101 est premier et a = 12 est primitif modulo 101.)
 14

Générateurs algorithmiques (pseudo-aléatoires, GPA)

Mini-exemple: On veut imiter des nombres de 1 à 100 tirés au hasard.

aft
1. Choisir un nombre x0 au hasard dans {1, . . . , 100}.
2. Pour n = 1, 2, 3, ..., retourner xn = 12 xn−1 mod 101 . Par exemple, si x0 = 1:

x1 = (12 × 1 mod 101) = 12,

x2 = (12 × 12 mod 101) = (144 mod 101) = 43,
x3 = (12 × 43 mod 101) = (516 mod 101) = 11, etc.
xn = 12n mod 101.

Visite tous les nombres de 1 à 100 une fois chacun avant de revenir à x0 .

u1
u2
u3
Dr
(Parce que m = 101 est premier et a = 12 est primitif modulo 101.)
Si on veut des nombres réels entre 0 et 1:
= x1 /101
= x2 /101
= x3 /101
=
=
=
12/101 ≈
43/101 ≈
11/101 ≈
0.11881188...,
0.42574257...,
0.10891089..., etc.

Ici, ui visite toutes les valeurs {1/101, 2/101, . . . , 100/101} quand i va de 1 à 100.
On a donc une excellente uniformité en une dimension.
 15

xn = 12 xn−1 mod 101; un = xn /101 = 12 un−1 mod 1

Paires de valeurs successives (un−1 , un ):

aft
(1/101, 12/101), (2/101, 24/101), (3/101, 36/101), ...
1

un

Dr 0 1 un−1
 15

xn = 12 xn−1 mod 101; un = xn /101 = 12 un−1 mod 1

Paires de valeurs successives (un−1 , un ):

aft
(1/101, 12/101), (2/101, 24/101), (3/101, 36/101), ...
1

un

Dr 0 1 un−1
 15

xn = 12 xn−1 mod 101; un = xn /101 = 12 un−1 mod 1

Paires de valeurs successives (un−1 , un ):

aft
(1/101, 12/101), (2/101, 24/101), (3/101, 36/101), ...
1

un

Dr 0 1 un−1
 15

xn = 12 xn−1 mod 101; un = xn /101 = 12 un−1 mod 1

Paires de valeurs successives (un−1 , un ):

aft
(1/101, 12/101), (2/101, 24/101), (3/101, 36/101), ...
1

un

Dr 0 1 un−1
 15

xn = 12 xn−1 mod 101; un = xn /101 = 12 un−1 mod 1

Paires de valeurs successives (un−1 , un ):

aft
(1/101, 12/101), (2/101, 24/101), (3/101, 36/101), ...
1

un

Dr 0 1 un−1
 16

xn = 51 xn−1 mod 101; un = xn /101.

1

aft
un

Dr 0 1 un−1

Ici, on a une bonne uniformité en une dimension, mais pas en deux!

On a les points (1/101, 51/101), (2/101, 1/101), (3/101, 52/101), ...
 17

Une récurrence linéaire plus longue

aft
On choisit 3 entiers x−2 , x−1 , x0 dans {0, 1, . . . , 4294967086}, pas tous 0. Puis pour
n = 1, 2, . . . , soit

xn = (1403580xn−2 − 810728xn−3 ) mod 4294967087,

un = xn /4294967087.

Dr
 17

Une récurrence linéaire plus longue

aft
On choisit 3 entiers x−2 , x−1 , x0 dans {0, 1, . . . , 4294967086}, pas tous 0. Puis pour
n = 1, 2, . . . , soit

xn = (1403580xn−2 − 810728xn−3 ) mod 4294967087,

un = xn /4294967087.

On peut prouver que la suite x0 , x1 , x2 , . . . est périodique, de période 42949670873 − 1 ≈ 296 ,

Dr
et que le vecteur (xn−2 , xn−1 , xn ) visite chacun des 42949670873 − 1 triplets non nuls
exactement une fois lorsque n parcourt un cycle.
 18

Combinaison de deux récurrences: MRG32k3a

On choisit 6 entiers:
x0 , x1 , x2 dans {0, 1, . . . , 4294967086} (pas tous 0) et

aft
y0 , y1 , y2 dans {0, 1, . . . , 4294944442} (pas tous 0).

xn = (1403580xn−2 − 810728xn−3 ) mod 4294967087,

yn = (527612yn−1 − 1370589yn−3 ) mod 4294944443,
un = [(xn − yn ) mod 4294967087]/4294967087.

Dr
 18

Combinaison de deux récurrences: MRG32k3a

On choisit 6 entiers:
x0 , x1 , x2 dans {0, 1, . . . , 4294967086} (pas tous 0) et

aft
y0 , y1 , y2 dans {0, 1, . . . , 4294944442} (pas tous 0).

xn = (1403580xn−2 − 810728xn−3 ) mod 4294967087,

yn = (527612yn−1 − 1370589yn−3 ) mod 4294944443,
un = [(xn − yn ) mod 4294967087]/4294967087.

(xn−2 , xn−1 , xn ) visite chacune des 42949670873 − 1 valeurs possibles.

(yn−2 , yn−1 , yn ) visite chacune des 42949444433 − 1 valeurs possibles.

Dr
La suite u0 , u1 , u2 , . . . se répète avec une période proche de 2191 ≈ 3.1 × 1057 .
L’uniformité des points (un+1 , . . . , un+s ) dans [0, 1)s a été mesurée mathématiquement
jusqu’en s = 48 dimensions.
 18

Combinaison de deux récurrences: MRG32k3a

On choisit 6 entiers:
x0 , x1 , x2 dans {0, 1, . . . , 4294967086} (pas tous 0) et

aft
y0 , y1 , y2 dans {0, 1, . . . , 4294944442} (pas tous 0).

xn = (1403580xn−2 − 810728xn−3 ) mod 4294967087,

yn = (527612yn−1 − 1370589yn−3 ) mod 4294944443,
un = [(xn − yn ) mod 4294967087]/4294967087.

(xn−2 , xn−1 , xn ) visite chacune des 42949670873 − 1 valeurs possibles.

(yn−2 , yn−1 , yn ) visite chacune des 42949444433 − 1 valeurs possibles.

Dr
La suite u0 , u1 , u2 , . . . se répète avec une période proche de 2191 ≈ 3.1 × 1057 .
L’uniformité des points (un+1 , . . . , un+s ) dans [0, 1)s a été mesurée mathématiquement
jusqu’en s = 48 dimensions.
Excellent générateur, robuste et fiable!
Disponible dans SSJ, SAS, R, MATLAB, Arena, Automod, Witness, machines Spielo, ...
 19

Générateurs algorithmiques
C’est ce qu’on utilise pour la simulation.

aft
Une fois les paramètres et l’état initial choisis, la suite devient complètement déterministe.
On peut choisir l’état initial au hasard si on veut.

Dr
 19

Générateurs algorithmiques
C’est ce qu’on utilise pour la simulation.

aft
Une fois les paramètres et l’état initial choisis, la suite devient complètement déterministe.
On peut choisir l’état initial au hasard si on veut.

Avantages: pas de matériel à installer, un logiciel suffit; souvent plus rapide; on peut
facilement répéter la même séquence.

Dr
 19

Générateurs algorithmiques
C’est ce qu’on utilise pour la simulation.

aft
Une fois les paramètres et l’état initial choisis, la suite devient complètement déterministe.
On peut choisir l’état initial au hasard si on veut.

Avantages: pas de matériel à installer, un logiciel suffit; souvent plus rapide; on peut
facilement répéter la même séquence.

Dr
Désavantage: ne peut pas créer de l’entropie!
Il y a nécessairement des dépendances entre les nombres en sortie.

Qualités requises pour un bon générateur: dépend des applications.

Pour la cryptographie, un générateur algorithmique seul ne suffit pas.
 20

Générateur algorithmique
S, espace d’états fini; s0 , germe (état initial);

aft
f : S → S, fonction de transition;
g : S → [0, 1], fonction de sortie.

s0

Dr
 20

Générateur algorithmique
S, espace d’états fini; s0 , germe (état initial);

aft
f : S → S, fonction de transition;
g : S → [0, 1], fonction de sortie.

s0

gy


u0

Dr
 20

Générateur algorithmique
S, espace d’états fini; s0 , germe (état initial);

aft
f : S → S, fonction de transition;
g : S → [0, 1], fonction de sortie.

f
s0 −−−−→ s1

gy


u0

Dr
 20

Générateur algorithmique
S, espace d’états fini; s0 , germe (état initial);

aft
f : S → S, fonction de transition;
g : S → [0, 1], fonction de sortie.

f
s0 −−−−→ s1
 
gy gy
 

u0 u1

Dr
 20

Générateur algorithmique
S, espace d’états fini; s0 , germe (état initial);

aft
f : S → S, fonction de transition;
g : S → [0, 1], fonction de sortie.

f f f f f
s0 −−−−→ s1 −−−−→ · · · −−−−→ sn −−−−→ sn+1 −−−−→ · · ·
   
gy gy gy gy
   

u0 u1 ··· un un+1 ···

Dr
 20

Générateur algorithmique
S, espace d’états fini; s0 , germe (état initial);

aft
f : S → S, fonction de transition;
g : S → [0, 1], fonction de sortie.

f f f f f f f
· · · −−−−→ sρ−1 −−−−→ s0 −−−−→ s1 −−−−→ · · · −−−−→ sn −−−−→ sn+1 −−−−→ · · ·
    
gy gy gy gy gy
    

··· uρ−1 u0 u1 ··· un un+1 ···

Dr
Période de {sn , n ≥ 0}: ρ ≤ cardinalité de S.

Exemple: Dans MRG32k3a, sn est un vecteur de 6 entiers de 32 bits et ρ ≈ 2191 .

 21

f f f f f f f
· · · −−−−→ sρ−1 −−−−→ s0 −−−−→ s1 −−−−→ · · · −−−−→ sn −−−−→ sn+1 −−−−→ · · ·
    

aft
gy gy gy gy gy
    

··· uρ−1 u0 u1 ··· un un+1 ···

Objectif: en observant seulement (u0 , u1 , . . .), difficile de distinguer d’une suite de v.a.
indépendantes uniformes sur (0, 1).

Dr
 21

f f f f f f f
· · · −−−−→ sρ−1 −−−−→ s0 −−−−→ s1 −−−−→ · · · −−−−→ sn −−−−→ sn+1 −−−−→ · · ·
    

aft
gy gy gy gy gy
    

··· uρ−1 u0 u1 ··· un un+1 ···

Objectif: en observant seulement (u0 , u1 , . . .), difficile de distinguer d’une suite de v.a.
indépendantes uniformes sur (0, 1).
Utopie: passe tous les tests statistiques imaginables.
Impossible! On doit se contenter d’une approximation.

Dr
 21

f f f f f f f
· · · −−−−→ sρ−1 −−−−→ s0 −−−−→ s1 −−−−→ · · · −−−−→ sn −−−−→ sn+1 −−−−→ · · ·
    

aft
gy gy gy gy gy
    

··· uρ−1 u0 u1 ··· un un+1 ···

Objectif: en observant seulement (u0 , u1 , . . .), difficile de distinguer d’une suite de v.a.
indépendantes uniformes sur (0, 1).
Utopie: passe tous les tests statistiques imaginables.
Impossible! On doit se contenter d’une approximation.

Dr
On veut aussi: vitesse, facilité d’implantation, suites reproduisibles.
Compromis entre vitesse / propriétés statistiques / imprévisibilité.
 21

f f f f f f f
· · · −−−−→ sρ−1 −−−−→ s0 −−−−→ s1 −−−−→ · · · −−−−→ sn −−−−→ sn+1 −−−−→ · · ·
    

aft
gy gy gy gy gy
    

··· uρ−1 u0 u1 ··· un un+1 ···

Objectif: en observant seulement (u0 , u1 , . . .), difficile de distinguer d’une suite de v.a.
indépendantes uniformes sur (0, 1).
Utopie: passe tous les tests statistiques imaginables.
Impossible! On doit se contenter d’une approximation.

Dr
On veut aussi: vitesse, facilité d’implantation, suites reproduisibles.
Compromis entre vitesse / propriétés statistiques / imprévisibilité.

Machines de casinos et loteries: on modifie l’état sn régulièrement à l’aide de mécanismes

physiques.
 22

La loi uniforme sur [0, 1]s .

aft
Si on choisit s0 au hasard dans S et on génère s nombres, cela correspond à choisir un point
au hasard dans l’ensemble fini

Ψs = {u = (u0 , . . . , us−1 ) = (g (s0 ), . . . , g (ss−1 )), s0 ∈ S}.

On veut approximer: “u suit la loi uniforme sur [0, 1]s .”

Dr
 22

La loi uniforme sur [0, 1]s .

aft
Si on choisit s0 au hasard dans S et on génère s nombres, cela correspond à choisir un point
au hasard dans l’ensemble fini

Ψs = {u = (u0 , . . . , us−1 ) = (g (s0 ), . . . , g (ss−1 )), s0 ∈ S}.

On veut approximer: “u suit la loi uniforme sur [0, 1]s .”

Mesure de qualité: Ψs doit recouvrir [0, 1]s très uniformément.

Dr
 22

La loi uniforme sur [0, 1]s .

aft
Si on choisit s0 au hasard dans S et on génère s nombres, cela correspond à choisir un point
au hasard dans l’ensemble fini

Ψs = {u = (u0 , . . . , us−1 ) = (g (s0 ), . . . , g (ss−1 )), s0 ∈ S}.

On veut approximer: “u suit la loi uniforme sur [0, 1]s .”

Mesure de qualité: Ψs doit recouvrir [0, 1]s très uniformément.

Dr
Conception et analyse théorique des générateurs:
1. Définir une mesure d’uniformité de Ψs , calculable
sans générer les points explicitement. GPA linéaires.
2. Choisir un type de construction (rapide, longue période, etc.)
et chercher des paramètres qui “optimisent” cette mesure.
 23

Mythe 1. Après au moins 60 ans à étudier les GPA et des milliers d’articles publiés, ce

aft
problème est certainement réglé et les GPA disponibles dans les logiciels populaires sont
certainement fiables.

Dr
 23

Mythe 1. Après au moins 60 ans à étudier les GPA et des milliers d’articles publiés, ce

aft
problème est certainement réglé et les GPA disponibles dans les logiciels populaires sont
certainement fiables.
Non.

Mythe 2. Dans votre logiciel favori, le générateur a une période supérieure à 21000 . Il est
donc certainement excellent!

Dr
 23

Mythe 1. Après au moins 60 ans à étudier les GPA et des milliers d’articles publiés, ce

aft
problème est certainement réglé et les GPA disponibles dans les logiciels populaires sont
certainement fiables.
Non.

Mythe 2. Dans votre logiciel favori, le générateur a une période supérieure à 21000 . Il est
donc certainement excellent!
Non.

Dr
Exemple 1. un = (n/21000 ) mod 1 pour n = 0, 1, 2, ....
Exemple 2. Lagged-Fibonacci, subtract-with-borrow.
 24

Suites (“streams”) et sous-suites multiples

Un seul GPA (monolithique) ne suffit pas. On a souvent besoin de plusieurs flux (ou

aft
suites, ou “streams”) “indépendants” de nombres aléatoires. Exemples:
▶ exécuter une simulation sur plusieurs processeurs en parallèle,
▶ Comparaison de systèmes avec valeurs aléatoires communes
(important pour analyse de sensibilité, estimation de dérivées, optimisation, ...).

Dr
 24

Suites (“streams”) et sous-suites multiples

Un seul GPA (monolithique) ne suffit pas. On a souvent besoin de plusieurs flux (ou

aft
suites, ou “streams”) “indépendants” de nombres aléatoires. Exemples:
▶ exécuter une simulation sur plusieurs processeurs en parallèle,
▶ Comparaison de systèmes avec valeurs aléatoires communes
(important pour analyse de sensibilité, estimation de dérivées, optimisation, ...).
Un logiciel développé au DIRO fournit de tels objets appelés RandomStream’s, qui agisssent
comme des GPA virtuels. On peut en créer autant qu’on veut.
On partitionne la suite en segments disjoints (streams) de longueur ν:

s0
stream 1
sν
Dr
stream 2
s2ν
stream 3
s3ν
stream 4
s4ν
stream 5
s5ν

Par exemple, pour MRG32k3a, la période est ρ ≈ 2191 et on a pris des segments de longueur
ν = 2127 . On peut donc avoir jusqu’à 264 RandomStream’s disjoints.
...
 25

Suites multiples en Java

RandomStream stream1 = new MRG32k3a();

aft
RandomStream stream2 = new MRG32k3a();
double u = [Link](); ....
[Link](); ....
[Link]();
[Link]();

Current
state
⇓
start
stream
..
Dr ..
start
substream
.
next
substream
.. .
 26

Sauter en avant
xn = (a1 xn−1 + · · · + ak xn−k ) mod m, un = xn /m.

aft
État à l’étape n: sn = xn = (xn−k+1 , . . . , xn )t .
 
0 1 ··· 0
 .. .. .. 
xn = Axn−1 mod m =  . . . xn−1 mod m.

0 0 ··· 1
ak ak−1 · · · a1

Dr
 26

Sauter en avant
xn = (a1 xn−1 + · · · + ak xn−k ) mod m, un = xn /m.

aft
État à l’étape n: sn = xn = (xn−k+1 , . . . , xn )t .
 
0 1 ··· 0
 .. .. .. 
xn = Axn−1 mod m =  . . . xn−1 mod m.

0 0 ··· 1
ak ak−1 · · · a1

xn+ν = Aν xn mod m = (Aν mod m)xn mod m.

Aν mod m =
(
Dr
On peut précalculer Aν mod m via

(Aν/2 mod m)(Aν/2 mod m) mod m

A(Aν−1 mod m) mod m

Fonctionne pour tout m > 1. Utilisé pour les grands m et aussi m = 2.

Rendu ici, 1 nov. 2024
si ν est pair;
si ν est impair.
Dr
aft
27
Dr
aft
28
 29

Variables aléatoires non uniformes: Inversion

Une v.a. X a la fonction de répartition (cdf) F si F (x) = P[X ≤ x] pour tout x ∈ R.

aft
Si U ∼ U(0, 1) (une v.a. uniforme sur (0, 1)) et

X = F −1 (U) = min{x | F (x) ≥ U},

alors X est une v.a. dont la cdf est F .

Preuve: P[X ≤ x] = P[F −1 (U) ≤ x] = P[U ≤ F (x)] = F (x).

U
Dr
Deux ingrédients nécessaires:
1. Un bon générateur uniforme pour générer U;
X
F (x)

2. Une formule ou une bonne approximation pour calculer rapidement F −1 (U).

 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.

Dr
 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.
Inversion:

U = 1 − exp[−(λX )α ]

Dr
 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.
Inversion:

U = 1 − exp[−(λX )α ]
exp[−(λX )α ] = 1 − U

Dr
 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.
Inversion:

U = 1 − exp[−(λX )α ]
exp[−(λX )α ] = 1 − U
(λX )α = − ln(1 − U)

Dr
 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.
Inversion:

U = 1 − exp[−(λX )α ]
exp[−(λX )α ] = 1 − U
(λX )α = − ln(1 − U)

Dr λX = [− ln(1 − U)]1/α
 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.
Inversion:

U = 1 − exp[−(λX )α ]
exp[−(λX )α ] = 1 − U
(λX )α = − ln(1 − U)

Dr λX
X
= [− ln(1 − U)]1/α
= [− ln(1 − U)]1/α /λ = F −1 (U).
 30

Exemple: Loi de Weibull

X ∼ Weibull(α, λ).

aft
F (x) = 1 − exp[−(λx)α ] pour x > 0.
Inversion:

U = 1 − exp[−(λX )α ]
exp[−(λX )α ] = 1 − U
(λX )α = − ln(1 − U)

Dr λX
X
= [− ln(1 − U)]1/α
= [− ln(1 − U)]1/α /λ = F −1 (U).

Cas particulier: si α = 1, X ∼ Exponentielle(λ).

Pour générer: X = F −1 (U) = − ln(1 − U)/λ.
 31

Exemple: Une loi discrète

aft
Soit P[X = i] = pi où p0 = 1/2, p2 = 3/8, p4 = 1/8, et pi = 0 ailleurs.
Inversion: retourner 0 si U < 1/2, 2 si 1/2 ≤ U < 7/8, et 4 si U ≥ 7/8.

1
7/8 F (x)
UU
1/2

Dr 0
0 2 4 x
x
 32

Exemple: Loi discrète

P[X = xi ] = pi pour i = 0, 1, . . . .
On a F (xi ) = p0 + · · · + pi et (disons) F (x−1 ) = 0.

aft
L’inversion retourne X = xi ssi F (xi ) ≥ U > F (xi−1 ).

1
F (xi ) F (x)
UU
F (xi−1 )

Dr 0
xi−1 xi x
x

Algorithme: Générer U ∼ U(0, 1), trouver I = min{i|F (xi ) ≥ U}, et retourner xI .

Mais trouver I en essayant i = 0, 1, 2, 3, . . . (recherche séquentielle) peut être très long.
 33

Recherche par index

On partitionne (0, 1) en c intervalles de longueur 1/c.

aft
On calcule et mémorise is = inf{i : F (xi ) ≥ s/c} pour s = 0, . . . , c − 1.
Si s = ⌊cU⌋, alors U ∈ [s/c, (s + 1)/c) et on a I := F −1 (U) ∈ {is , . . . , is+1 }.
On calcule et mémorise aussi les F (xi ) pour i = 0, 1, 2, . . . .

Dr
 33

Recherche par index

On partitionne (0, 1) en c intervalles de longueur 1/c.

aft
On calcule et mémorise is = inf{i : F (xi ) ≥ s/c} pour s = 0, . . . , c − 1.
Si s = ⌊cU⌋, alors U ∈ [s/c, (s + 1)/c) et on a I := F −1 (U) ∈ {is , . . . , is+1 }.
On calcule et mémorise aussi les F (xi ) pour i = 0, 1, 2, . . . .

Il suffit alors de chercher dans cet intervalle, par recherche séquentielle ou binaire.
Inversion avec recherche par index (combiné avec recherche séquentielle);
générer U ∼ U(0, 1); poser s = ⌊cU⌋ et i = is ;
tant que F (xi ) < U faire i = i + 1;
retourner xi .

Dr
 33

Recherche par index

On partitionne (0, 1) en c intervalles de longueur 1/c.

aft
On calcule et mémorise is = inf{i : F (xi ) ≥ s/c} pour s = 0, . . . , c − 1.
Si s = ⌊cU⌋, alors U ∈ [s/c, (s + 1)/c) et on a I := F −1 (U) ∈ {is , . . . , is+1 }.
On calcule et mémorise aussi les F (xi ) pour i = 0, 1, 2, . . . .

Il suffit alors de chercher dans cet intervalle, par recherche séquentielle ou binaire.
Inversion avec recherche par index (combiné avec recherche séquentielle);
générer U ∼ U(0, 1); poser s = ⌊cU⌋ et i = is ;
tant que F (xi ) < U faire i = i + 1;
retourner xi .

Dr
Le nombre espéré d’itérations du “tant que” est environ k/c.
Si on choisit c ≈ k, par exemple, alors on a un algorithme super rapide.
On paye un peu en terme de mémoire.
En pratique, on pourra prendre c jusqu’à 212 = 4096, par exemple.
 34

Recherche par index

Posons s = ⌊cU⌋.

aft
On a U ∈ [s/c, (s + 1)/c) et I := F −1 (U) ∈ {is , . . . , is+1 }.

1
F (x)
(s + 1)/c = 3/4
UU
s/c = 2/4

Dr 0
x4 x5 x6 x7

Pour ce U, on a c = 4, s = 2, is = 4, is+1 = 7, donc I ∈ {4, 5, 6, 7}.

x
x
 35

Exemple: Inversion pour la loi géométrique

X ∼ Géométrique(p), où 0 < p < 1.

aft
P[X = x] = p(1 − p)x pour x = 0, 1, 2, . . . et
F (x) = 1 − (1 − p)⌊x+1⌋ pour x ≥ 0.

Dr
 35

Exemple: Inversion pour la loi géométrique

X ∼ Géométrique(p), où 0 < p < 1.

aft
P[X = x] = p(1 − p)x pour x = 0, 1, 2, . . . et
F (x) = 1 − (1 − p)⌊x+1⌋ pour x ≥ 0.

Pour x ≥ 0 entier, F (x) = 1 − (1 − p)x+1 .

L’inversion doit retourner X = x ssi F (x) ≥ U > F (x − 1),

Dr
 35

Exemple: Inversion pour la loi géométrique

X ∼ Géométrique(p), où 0 < p < 1.

aft
P[X = x] = p(1 − p)x pour x = 0, 1, 2, . . . et
F (x) = 1 − (1 − p)⌊x+1⌋ pour x ≥ 0.

Pour x ≥ 0 entier, F (x) = 1 − (1 − p)x+1 .

L’inversion doit retourner X = x ssi F (x) ≥ U > F (x − 1),
que l’on peut réécrire

1 − (1 − p)x+1 ≥ U > 1 − (1 − p)x ,

Dr
 35

Exemple: Inversion pour la loi géométrique

X ∼ Géométrique(p), où 0 < p < 1.

aft
P[X = x] = p(1 − p)x pour x = 0, 1, 2, . . . et
F (x) = 1 − (1 − p)⌊x+1⌋ pour x ≥ 0.

Pour x ≥ 0 entier, F (x) = 1 − (1 − p)x+1 .

L’inversion doit retourner X = x ssi F (x) ≥ U > F (x − 1),
que l’on peut réécrire

1 − (1 − p)x+1 ≥ U > 1 − (1 − p)x , ou

Dr(1 − p)x+1 ≤ 1 − U < (1 − p)x ,

 35

Exemple: Inversion pour la loi géométrique

X ∼ Géométrique(p), où 0 < p < 1.

aft
P[X = x] = p(1 − p)x pour x = 0, 1, 2, . . . et
F (x) = 1 − (1 − p)⌊x+1⌋ pour x ≥ 0.

Pour x ≥ 0 entier, F (x) = 1 − (1 − p)x+1 .

L’inversion doit retourner X = x ssi F (x) ≥ U > F (x − 1),
que l’on peut réécrire

1 − (1 − p)x+1 ≥ U > 1 − (1 − p)x , ou

Dr(1 − p)x+1 ≤ 1 − U < (1 − p)x ,

(x + 1) ln(1 − p) ≤ ln(1 − U) < x ln(1 − p),
ou
 35

Exemple: Inversion pour la loi géométrique

X ∼ Géométrique(p), où 0 < p < 1.

aft
P[X = x] = p(1 − p)x pour x = 0, 1, 2, . . . et
F (x) = 1 − (1 − p)⌊x+1⌋ pour x ≥ 0.

Pour x ≥ 0 entier, F (x) = 1 − (1 − p)x+1 .

L’inversion doit retourner X = x ssi F (x) ≥ U > F (x − 1),
que l’on peut réécrire

1 − (1 − p)x+1 ≥ U > 1 − (1 − p)x , ou

Dr
(1 − p)x+1 ≤ 1 − U < (1 − p)x ,
(x + 1) ln(1 − p) ≤ ln(1 − U) < x ln(1 − p),
(x + 1) ≥ ln(1 − U)/ ln(1 − p) > x, (car ln(1 − p) < 0)

ce qui donne x + 1 = ⌈ln(1 − U)/ ln(1 − p)⌉.

ou
ou

Avec probabilité 1, c’est la même chose que retourner X = ⌊ln(1 − U)/ ln(1 − p)⌋.
 36

Autres situations
▶ Dans plusieurs cas (normale, Student, chi-deux, etc.), pas de formule pour F −1 , mais

aft
approximation numérique.

Dr
 36

Autres situations
▶ Dans plusieurs cas (normale, Student, chi-deux, etc.), pas de formule pour F −1 , mais

aft
approximation numérique.
▶ Plus difficile lorsque la forme de F dépend des paramètres (beta, gamma, par exemple).

Dr
 36

Autres situations
▶ Dans plusieurs cas (normale, Student, chi-deux, etc.), pas de formule pour F −1 , mais

aft
approximation numérique.
▶ Plus difficile lorsque la forme de F dépend des paramètres (beta, gamma, par exemple).
▶ Inversion préférable car monotone (on verra pourquoi plus tard).

Dr
 36

Autres situations
▶ Dans plusieurs cas (normale, Student, chi-deux, etc.), pas de formule pour F −1 , mais

aft
approximation numérique.
▶ Plus difficile lorsque la forme de F dépend des paramètres (beta, gamma, par exemple).
▶ Inversion préférable car monotone (on verra pourquoi plus tard).
▶ Mais d’autres méthodes sont parfois beaucoup plus rapides.

Dr
 36

Autres situations
▶ Dans plusieurs cas (normale, Student, chi-deux, etc.), pas de formule pour F −1 , mais

aft
approximation numérique.
▶ Plus difficile lorsque la forme de F dépend des paramètres (beta, gamma, par exemple).
▶ Inversion préférable car monotone (on verra pourquoi plus tard).
▶ Mais d’autres méthodes sont parfois beaucoup plus rapides.

Inversion pour la loi de Y conditionnelle à Y ∈ (a, b].

Théorème: Soit F la cdf d’une variable aléatoire Y et soit −∞ ≤ a < b ≤ ∞.
Si U ∼ Uniforme(F (a), F (b)) et X = F −1 (U), alors la loi de probabilité de X (sans

Dr
condition) est la même que la loi de Y conditionnelle à Y ∈ (a, b].
Preuve: Puisque F (X ) = U, on a
P(X ≤ x) = P(U ≤ F (x)) = FF (x)−F (a)
(b)−F (a) pour F (a) ≤ F (x) ≤ F (b).
C’est la même cdf que la cdf de Y conditionnelle à Y ∈ (a, b]:
P[Y ≤ y | a < Y ≤ b] = FF (y )−F (a)
(b)−F (a) pour a ≤ y ≤ b.
 37

Méthode de rejet
Technique la plus importante après l’inversion.
Peut fournir une solution efficace lorsque l’inversion est trop difficile ou coûteuse.

aft
On veut générer X selon une densité f . La surface sous f est:

S(f ) = {(x, y ) ∈ R2 : 0 ≤ y ≤ f (x)}.

Proposition. Si le point (X , Y ) est uniforme sur S(f ), alors X a la densité f .

R x S(f ), alors P[X ≤ x] est égal à la surface de

Preuve. Si (X , Y ) est uniforme sur
{(z, y ) ∈ S(f ) : z ≤ x}, qui est −∞ f (z)dz. □

Dr
 37

Méthode de rejet
Technique la plus importante après l’inversion.
Peut fournir une solution efficace lorsque l’inversion est trop difficile ou coûteuse.

aft
On veut générer X selon une densité f . La surface sous f est:

S(f ) = {(x, y ) ∈ R2 : 0 ≤ y ≤ f (x)}.

Proposition. Si le point (X , Y ) est uniforme sur S(f ), alors X a la densité f .

R x S(f ), alors P[X ≤ x] est égal à la surface de

Preuve. Si (X , Y ) est uniforme sur
{(z, y ) ∈ S(f ) : z ≤ x}, qui est −∞ f (z)dz. □

Dr
On va générer (X , Y ) uniformément sur S(f ). Comment, si S(f ) est compliqué?
 37

Méthode de rejet
Technique la plus importante après l’inversion.
Peut fournir une solution efficace lorsque l’inversion est trop difficile ou coûteuse.

aft
On veut générer X selon une densité f . La surface sous f est:

S(f ) = {(x, y ) ∈ R2 : 0 ≤ y ≤ f (x)}.

Proposition. Si le point (X , Y ) est uniforme sur S(f ), alors X a la densité f .

R x S(f ), alors P[X ≤ x] est égal à la surface de

Preuve. Si (X , Y ) est uniforme sur
{(z, y ) ∈ S(f ) : z ≤ x}, qui est −∞ f (z)dz. □

Dr
On va générer (X , Y ) uniformément sur S(f ). Comment, si S(f ) est compliqué?

Idée: Choisir une surface simple B qui contient S(f ), et générer (X , Y ) uniformément dans
B. Si (X , Y ) ∈ S(f ), c’est bon, sinon on recommence.
On va montrer que le point (X , Y ) retenu suit une loi uniforme sur S(f ).
 38
Exemple: On veut générer X ∼ Beta(3, 2), de densité f (x) = 12x 2 (1 − x) sur (0, 1).
La densité est maximale à x∗ = 2/3. On a f (2/3) = 16/9 ≈ 1.77778.
Ici, S(f ) est la surface sous la courbe en rouge.

aft
Alors on peut prendre B = {(x, y ) : 0 ≤ x ≤ 1, 0 ≤ y ≤ 16/9} (un rectangle).
h(x)
↓
a = 16/9 · · · · · · · · · · · · ········· · · · · · ·
....................................................................................................................................................................................................................................................................................................
....
....
..........
.........
· · · · · · · · ·· ································ · · ·
.......
... ........ ......
....... ......
... .......
· · · · · aV · · · ······························· · · · . ............
.....

· · · · · · →
. .....
• .....
· ················· · · ........
.
. .....
1 · · · · · · · · · · · · ················ ········································································ · · ·
..
......... .....
........ ...
. ...
...... .
. ...

· · · · · ······························ ························································· ↙· f (x)

... . ...
...
...... .
.
. ...
......... .
.
.
...
... ...
· · · · ·········· · · · · · · · · · · · · · · · · · · · · ...
. .
.
· · · · ························································ ····························································· ·
.
...... .
. ...
.
.... .
. ...
.......... .
.
. ...
..... ...
· · ················ · · · · · · · · · · · · · · · · · · · ·
...
. .
.
· · · ··················································· ···········································
......... .
. ...
...
.. .
0 ····· ..... . ...
...................

Dr
(X , Y ) = (U, aV ) où a = 16/9.
La probabilité que le point soit dans S(f ) est 1/a = 9/16.
Le nombre espéré de points (X , Y ) qu’il faudra générer est a = 16/9.
X

Pour générer un point dans B, on génère U, V ∼ U(0, 1) indépendantes, et on pose

2/3 1
.
 39

Méthode de rejet générale

aft
On veut générer un point X uniformément dans un ensemble A ⊂ Rd .
On choisit un ensemble B plus “simple” tel que A ⊂ B.
On génère des points indépendants dans B, et on retient le premier qui tombe dans A.
Proposition. Le point X retenu suit la loi uniforme sur A.
Preuve. On veut montrer que pour tout D ⊆ A, on a P[X ∈ D | X ∈ A] = vol(D)/vol(A).
Puisque X est uniforme sur B, on a P[X ∈ D] = vol(D)/vol(B) et donc

Dr
P[X ∈ D | X ∈ A] =
P[X ∈ D ∩ A]
P[X ∈ A]
=
vol(D)/vol(B)
vol(A)/vol(B)

Ainsi la loi de X conditionelle à X ∈ A est uniforme sur A.

=
vol(D)
vol(A)
.
 40

Méthode de rejet avec fonction chapeau

Pour générer X selon f , on choisit une autre densité g et une constante a ≥ 1 telle que

aft
def
f (x) ≤ h(x) = ag (x)

pour tout x, et telle qu’il est facile de générer X selon g . La fonction chapeau est h.

Dr
 40

Méthode de rejet avec fonction chapeau

Pour générer X selon f , on choisit une autre densité g et une constante a ≥ 1 telle que

aft
def
f (x) ≤ h(x) = ag (x)

pour tout x, et telle qu’il est facile de générer X selon g . La fonction chapeau est h.
On applique la méthode de rejet avec A = S(f ) et

B = S(h) = {(x, y ) ∈ R2 : 0 ≤ y ≤ h(x)}, la surface sous h.

Algorithme de rejet;
répéter

Dr
générer X selon la densité g et V ∼ U(0, 1), indépendants;
// Le point (X , h(X )V ) est uniforme sur S(h).
jusqu’à ce que V h(X ) ≤ f (X );
retourner X . // Le point X retourné est uniforme sur A = S(f ).
Proposition. La v.a. X retournée a la densité f .
 41

aft
Algorithme de rejet;
répéter
générer X selon la densité g et V ∼ U(0, 1), indépendants;
jusqu’à ce que V h(X ) ≤ f (X );
retourner X .

À chaque tour de boucle, la probabilité d’accepter X est 1/a.

Le nombre R de rejets avant l’acceptation est une v.a. géométrique de paramètre p = 1/a.

Dr
Le nombre moyen de tours de boucle requis par v.a. est donc 1/p = a.
On veut a ≥ 1 le plus petit possible.
Compromis entre diminuer a et garder g simple.
 42

Exemple: X ∼ Beta(3, 2) (suite). Baisser le chapeau.

Pour réduire a, on peut prendre la fonction chapeau:

aft

f (x1 ) pour x < x1 ;

h(x) = 16/9 pour x1 ≤ x ≤ x2 ;

f (x2 ) pour x > x2 ,


où x1 et x2 satisfont 0 < x1 < 2/3 < x2 < 1.

h(x)
↓
a = 16/9

0
Dr 0
R1

..............
·

...
· · · · ····················· · ·
...........................................................................................................................................................................................................................................................................................................................................................

...
· · ············· ·
· · ····························· · · R2
· ···························································
· · · · · · · · · · · · · · · · · · · · · · · · f (x)
· · · · · · ··············································· ↙
..... .
.......
·····································
. ..............................
......

· · ·······························································
.
..
..
...
...

·
· · · ·· ················································································ ·
...

.........................................................................................

x1
...
...
... ............
.
......
......
. .
..... ..... ..... ..........
.........
..........
.........

2/3
..........
........
.......
......
........ ....

x2
..
...

........
. .
...
...

................................
.....
.....
....
...
...

1
...
...
...
...
 43

aft
On génère un point au hasard dans la surface sous h et on l’accepte s’il est sous la courbe
rouge.
La surface sous h est minimisée en prenant x1 = 0.281023 et x2 = 0.89538. Elle est alors
réduite de 1.77778 à 1.38997. La probabilité d’accepter passe de 1/1.77778 à 1/1.38997.
La fonction de répartition inverse de g est linéaire par morceaux.
Pourquoi ne pas prendre une fonction h linéaire par morceaux au lieu de constante par

Dr
morceaux? Le calcul de G −1 demande alors des racines carrées... et on perd en efficacité.
 44

Point au hasard sur la surface d’une hypersphère

Sphère de rayon 1 centrée à l’origine, en d dimensions.

aft
Il suffit de généner un point dans Rd , selon une densité radialement symétrique (i.e., qui ne
dépend que de la distance à l’origine). Par exemple, la densité multinormale standard.

Dr
 44

Point au hasard sur la surface d’une hypersphère

Sphère de rayon 1 centrée à l’origine, en d dimensions.

aft
Il suffit de généner un point dans Rd , selon une densité radialement symétrique (i.e., qui ne
dépend que de la distance à l’origine). Par exemple, la densité multinormale standard.

Pd
On génère Z = (Z1 , . . . , Zd ) ∼ N(0, I ) et on pose X = Z /∥Z ∥2 , où ∥Z ∥22 = 2
j=1 Zj .
Ici, Z1 , . . . , Zd sont indépendantes et N(0, 1).

Dr
 44

Point au hasard sur la surface d’une hypersphère

Sphère de rayon 1 centrée à l’origine, en d dimensions.

aft
Il suffit de généner un point dans Rd , selon une densité radialement symétrique (i.e., qui ne
dépend que de la distance à l’origine). Par exemple, la densité multinormale standard.

Pd
On génère Z = (Z1 , . . . , Zd ) ∼ N(0, I ) et on pose X = Z /∥Z ∥2 , où ∥Z ∥22 = 2
j=1 Zj .
Ici, Z1 , . . . , Zd sont indépendantes et N(0, 1).

En deux dimensions: on peut générer Θ uniformément sur (0, 2π), puis poser
(X1 , X2 ) = (cos Θ, sin Θ).

Dr
 44

Point au hasard sur la surface d’une hypersphère

Sphère de rayon 1 centrée à l’origine, en d dimensions.

aft
Il suffit de généner un point dans Rd , selon une densité radialement symétrique (i.e., qui ne
dépend que de la distance à l’origine). Par exemple, la densité multinormale standard.

Pd
On génère Z = (Z1 , . . . , Zd ) ∼ N(0, I ) et on pose X = Z /∥Z ∥2 , où ∥Z ∥22 = 2
j=1 Zj .
Ici, Z1 , . . . , Zd sont indépendantes et N(0, 1).

En deux dimensions: on peut générer Θ uniformément sur (0, 2π), puis poser
(X1 , X2 ) = (cos Θ, sin Θ).

r˜ = cos(arcsin x3 ).
Dr
En trois dimensions: chaque coordonnée Xj est uniforme sur (−1, 1) !
On peut générer X3 ∼ U(−1, 1), puis (X1 , X2 ) sur le cercle qui reste, dont le rayon est

Pour cela, générer Θ ∼ U(0, 2π), et poser (X1 , X2 ) = (˜

r cos Θ, r˜ sin Θ).
 45

Méthode Monte Carlo

Intégration par Monte Carlo: Monte Carlo est souvent utilisé pour estimer une espérance

aft
µ = E[X ]

où X est une v.a. dont chaque réalisation est calculée par simulation. On peut interpréter X
comme une fonction d’un vecteur U = (U0 , . . . , Us−1 ) de s v.a. i.i.d. U(0, 1),
X = f (U) = f (U1 , . . . , Us ), et on a
Z Z 1 Z 1
µ = E[X ] = E[f (U)] = f (u)du = ··· f (u0 , . . . , us−1 )du0 · · · dus−1 ,

Dr [0,1)s 0

La dimension s peut parfois être aléatoire ou infinie.

La fonction f est souvent très compliquée.
µ est la valeur moyenne de f sur [0, 1)s .
Rendu ici, 4 nov. 2024
0
 46

aft
Exemple: réseau d’activités stochastique à 13 activités. On a s = 13.
À un point U = (U0 , . . . , U12 ) ∈ [0, 1)13 correspond un vecteur de durées (Y0 , . . . , Y12 ) =
(F0−1 (U0 ), . . . , F12
−1
(U12 )) et une durée du projet T .

Dr
 46

aft
Exemple: réseau d’activités stochastique à 13 activités. On a s = 13.
À un point U = (U0 , . . . , U12 ) ∈ [0, 1)13 correspond un vecteur de durées (Y0 , . . . , Y12 ) =
(F0−1 (U0 ), . . . , F12
−1
(U12 )) et une durée du projet T .

Pour estimer µ = E[T ], on peut définir f (U0 , . . . , U12 ) = T .

Pour estimer µ = P[T > x] = E[I(T > x)], on peut définir f (U0 , . . . , U12 ) = I(T > x).
Rappel: I(T > x) = 1 si T > x et = 0 sinon.

Dr
Dans les deux cas, on a une intégrale à 13 dimensions.
 47
Si s ne dépasse pas 2 ou 3 et f est assez lisse, peut utiliser des méthodes d’intégration
numériques (méthode du trapèze, méthode de Simpson, etc.).

aft
Dr
 47
Si s ne dépasse pas 2 ou 3 et f est assez lisse, peut utiliser des méthodes d’intégration
numériques (méthode du trapèze, méthode de Simpson, etc.).

aft
Si s est grand, estimateur de Monte Carlo de µ = E[X ] pour X = f (U):
n
1X
µ̂n = Xi où Xi = f (Ui ),
n
i=1

n est la taille de l’échantillon (nombre de simulations) et U1 , . . . , Un sont i.i.d. uniformes sur

[0, 1)s .
Z

Dr
E[µ̂n ] = E[f (Ui )] =

Var[µ̂n ] =

où σ 2 =
1 2
Zn
f (u)du = µ = E[X ]
[0,1)s
 1
E[f (Ui )] − µ2 = Var[X ] =

[0,1)s
n
σ2
n
f 2 (u)du − µ2 = Var[f (Ui )] = Var[X ]
 48

Convergence
Théorème. Supposons que σ 2 < ∞. Lorsque n → ∞:
(i) Loi forte des grands nombres: limn→∞ µ̂n = µ avec probabilité 1.

aft
Dr
 48

Convergence
Théorème. Supposons que σ 2 < ∞. Lorsque n → ∞:
(i) Loi forte des grands nombres: limn→∞ µ̂n = µ avec probabilité 1.

aft
(ii) Théorème de la limite centrale:
√
n(µ̂n − µ)
⇒ N(0, 1), i.e.,
σ √
n(µ̂n − µ)

lim P ≤ x = Φ(x) = P[Z ≤ x]
n→∞ σ

pour tout x ∈ R, où Z ∼ N(0, 1) et Φ(·) sa fonction de répartition.

Dr
 48

Convergence
Théorème. Supposons que σ 2 < ∞. Lorsque n → ∞:
(i) Loi forte des grands nombres: limn→∞ µ̂n = µ avec probabilité 1.

aft
(ii) Théorème de la limite centrale:
√
n(µ̂n − µ)
⇒ N(0, 1), i.e.,
σ √
n(µ̂n − µ)

lim P ≤ x = Φ(x) = P[Z ≤ x]
n→∞ σ

pour tout x ∈ R, où Z ∼ N(0, 1) et Φ(·) sa fonction de répartition.

2
Sn =

où Xi = f (Ui ). On a
Dr
La propriété (ii) tient aussi si on remplace σ 2 par son estimateur sans biais

√
n−1
1 X
n

i=1
2
(Xi − X̄n ) =
1
n−1
X n

i=1
2 2
Xi − n(X̄n ) ,
!

n(µ̂n − µ)/Sn ⇒ N(0, 1) . Donne une idée de la distrib. de l’erreur.

 49
Pour n grand et un niveau de confiance 1 − α, on a
√ √
P[µ̂n − µ ≤ xSn / n] = P[ n(µ̂n − µ)/Sn ≤ x] ≈ Φ(x).

aft
Intervalle de confiance au niveau 1 − α (on veut Φ(x) = 1 − α/2):
√
(µ̂n ± zα/2 Sn / n), où Φ(zα/2 ) = 1 − α/2.
Exemple: zα/2 ≈ 1.96 pour α = 0.05.

−zα/2 zα/2
α/2 1−α α/2

−3 −1.96 −1 0 1 1.96 3

Dr √
La largeur de l’intervalle de confiance est asymptotiquement proportionnelle à σ/ n, donc
converge en O(n−1/2 ).
 49
Pour n grand et un niveau de confiance 1 − α, on a
√ √
P[µ̂n − µ ≤ xSn / n] = P[ n(µ̂n − µ)/Sn ≤ x] ≈ Φ(x).

aft
Intervalle de confiance au niveau 1 − α (on veut Φ(x) = 1 − α/2):
√
(µ̂n ± zα/2 Sn / n), où Φ(zα/2 ) = 1 − α/2.
Exemple: zα/2 ≈ 1.96 pour α = 0.05.

−zα/2 zα/2
α/2 1−α α/2

−3 −1.96 −1 0 1 1.96 3

Dr √
√
La largeur de l’intervalle de confiance est asymptotiquement proportionnelle à σ/ n, donc
converge en O(n−1/2 ).
Si n est petit et les Xi suivent la loi normale, alors n(µ̂n − µ)/Sn ∼ Student(n − 1).
Si la loi des Xi est très asymétrique, ces intervalles ne sont plus valides même comme
approximations.
 50

Exemple: Réseau d’activités stochastique. Pour µ = P[T > x], on a

aft
Xi = I[Ti > x],
n
1X Y
µ̂n = X̄n = Xi = ,
n n
i=1
n
1 X Y (1 − Y /n)
Sn2 = (Xi − Y /n)2 = .
n−1 n−1
i=1

Dr √
Intervalle de confiance à 95%: (Y /n ± 1.96Sn / n).
Raisonnable si µ n’est pas trop proche de 0 ou 1.
En fait, Y = ni=1 Xi ∼ Binomiale(n, µ).
P
 51

Supposons par ex. que n = 1000 et qu’on observe Y = 882. On a alors:

aft
X̄n = 882/1000 = 0.882;
Sn2 = X̄n (1 − X̄n )n/(n − 1) ≈ 0.1042.

Dr
 51

Supposons par ex. que n = 1000 et qu’on observe Y = 882. On a alors:

aft
X̄n = 882/1000 = 0.882;
Sn2 = X̄n (1 − X̄n )n/(n − 1) ≈ 0.1042.
On obtient l’intervalle de confiance à 95%:
√
(X̄n ± 1.96Sn / n) ≈ (0.882 ± 0.020) = (0.862, 0.902).

Dr
 51

Supposons par ex. que n = 1000 et qu’on observe Y = 882. On a alors:

aft
X̄n = 882/1000 = 0.882;
Sn2 = X̄n (1 − X̄n )n/(n − 1) ≈ 0.1042.
On obtient l’intervalle de confiance à 95%:
√
(X̄n ± 1.96Sn / n) ≈ (0.882 ± 0.020) = (0.862, 0.902).

Notre estimateur de µ a donc deux chiffres significatifs: µ ≈ 0.88.

Le “2” dans 0.882 n’est pas significatif.

Dr
 51

Supposons par ex. que n = 1000 et qu’on observe Y = 882. On a alors:

aft
X̄n = 882/1000 = 0.882;
Sn2 = X̄n (1 − X̄n )n/(n − 1) ≈ 0.1042.
On obtient l’intervalle de confiance à 95%:
√
(X̄n ± 1.96Sn / n) ≈ (0.882 ± 0.020) = (0.862, 0.902).

Notre estimateur de µ a donc deux chiffres significatifs: µ ≈ 0.88.

Le “2” dans 0.882 n’est pas significatif.

Dr
Il faut éviter de donner des chiffres non significatifs dans les résultats, car cela peut induire
en erreur.
 52

Avantages de MC pour l’intégration

aft
▶ Ne requiert qu’une hypothèse très faible sur f .
▶ Le taux de convergence de l’erreur ne dépend pas de la dimension s, contrairement aux
méthodes d’intégration numérique classiques.
▶ On peut estimer l’erreur de manière probabiliste (intervalles de confiance).
Les méthodes d’intégration numériques donnent des bornes déterministes sur l’erreur,

Dr
mais souvent on ne peut pas les calculer, donc peu pratiques.
 53

Efficacité des estimateurs en simulation

aft
Soit X estimateur de µ.
La variance Var[X ] n’est pas la seule mesure de qualité de l’estimateur X .

β = E[X ] − µ biais
2
σ = Var(X ) = E[(X − E[X ])2 ] variance
2 2 2
MSE[X ] = E[(X − µ) ] = β + σ erreur quadratique moyenne

RE[X ] =
p
p
Dr
MSE[X ] erreur absolue
MSE[X ]/|µ|, pour µ ̸= 0 erreur relative
 54

Soit C (X ) l’espérance mathématique du temps de calcul de X .

Variance normalisée par le coût: C (X ) · MSE(X ).

aft
Efficacité de l’estimateur X :
1
Eff(X ) = .
C (X ) · MSE(X )

X est plus efficace que Y si Eff(X ) > Eff(Y ).

Dr
 54

Soit C (X ) l’espérance mathématique du temps de calcul de X .

Variance normalisée par le coût: C (X ) · MSE(X ).

aft
Efficacité de l’estimateur X :
1
Eff(X ) = .
C (X ) · MSE(X )

X est plus efficace que Y si Eff(X ) > Eff(Y ).

Amélioration de l’efficacité: trouver des estimateurs plus efficaces, en ce sens, soit en
diminuant la variance, ou le biais, ou le temps de calcul.

Dr
 54

Soit C (X ) l’espérance mathématique du temps de calcul de X .

Variance normalisée par le coût: C (X ) · MSE(X ).

aft
Efficacité de l’estimateur X :
1
Eff(X ) = .
C (X ) · MSE(X )

X est plus efficace que Y si Eff(X ) > Eff(Y ).

Amélioration de l’efficacité: trouver des estimateurs plus efficaces, en ce sens, soit en
diminuant la variance, ou le biais, ou le temps de calcul.

Dr
Soit X1 , . . . , Xn i.i.d., E[Xi ] = µ, et C (X̄n ) = κn. On a Var[X̄n ] = Var[Xi ]/n = σ 2 /n et

Eff[X̄n ] =
1
C (X̄n ) · MSE(X̄n )
=
1
κnσ 2 /n

Cette mesure d’efficacité ne dépend pas de n, ce qui est bien.

=
1
κσ 2
.
 55

Exemple de difficulté: événements rares

On veut estimer p = P{A} où A est un événement rare (p est proche de 0).

aft
La variable binaire X = I[A] est un estimateur sans biais de p,
de variance (et MSE) Var[X ] = MSE[X ] = p(1 − p).

Dr
 55

Exemple de difficulté: événements rares

On veut estimer p = P{A} où A est un événement rare (p est proche de 0).

aft
La variable binaire X = I[A] est un estimateur sans biais de p,
de variance (et MSE) Var[X ] = MSE[X ] = p(1 − p).
Si p est petit, Var[X ] = MSE[X ] ≈ p est petite.

Dr
 55

Exemple de difficulté: événements rares

On veut estimer p = P{A} où A est un événement rare (p est proche de 0).

aft
La variable binaire X = I[A] est un estimateur sans biais de p,
de variance (et MSE) Var[X ] = MSE[X ] = p(1 − p).
Si p est petit, Var[X ] = MSE[X ] ≈ p est petite.
Mais l’estimateur trivial Y = 0 donne déjà Var[Y ] = 0 et MSE[Y ] = p 2 ,
ce qui est encore plus petit et ne coûte rien! On veut faire mieux que cela.

Si on prend X̄n comme estimateur, alors Var[X̄n ] = p(1 − p)/n.

On a MSE[X̄n ] < MSE[Y ] ssi p(1 − p)/n < p 2 ssi n > (1 − p)/p.

Dr
 55

Exemple de difficulté: événements rares

On veut estimer p = P{A} où A est un événement rare (p est proche de 0).

aft
La variable binaire X = I[A] est un estimateur sans biais de p,
de variance (et MSE) Var[X ] = MSE[X ] = p(1 − p).
Si p est petit, Var[X ] = MSE[X ] ≈ p est petite.
Mais l’estimateur trivial Y = 0 donne déjà Var[Y ] = 0 et MSE[Y ] = p 2 ,
ce qui est encore plus petit et ne coûte rien! On veut faire mieux que cela.

Si on prend X̄n comme estimateur, alors Var[X̄n ] = p(1 − p)/n.

On a MSE[X̄n ] < MSE[Y ] ssi p(1 − p)/n < p 2 ssi n > (1 − p)/p.

Dr
Si |p| est petit, il est plus approprié de considérer le MSE relatif MSE[X ]/p 2 , ou l’erreur
relative RE[X ], car la largeur relative d’un intervalle depconfiance sur p est à peu près
√
proportionnelle à RE[X̄n ] = RE[X ]/ n. Ici, RE[X ] = (1 − p)/p → ∞ lorsque p → 0.
Par ex., si p ≈ 10−10 , il faut n ≈ 1012 pour avoir RE[X̄n ] ≈ 10%.
 56

aft
Autre point de vue: Si p est très petit (np est petit), on a de fortes chances d’avoir
X1 = · · · = Xn = 0, ce qui donne X̄n = Sn2 = 0 et un intervalle de confiance de largeur 0.

P[X̄n = Sn2 = 0] = (P[Xi = 0])n = (1 − p)n ≈ 1 − np.

Dr
 56

aft
Autre point de vue: Si p est très petit (np est petit), on a de fortes chances d’avoir
X1 = · · · = Xn = 0, ce qui donne X̄n = Sn2 = 0 et un intervalle de confiance de largeur 0.

P[X̄n = Sn2 = 0] = (P[Xi = 0])n = (1 − p)n ≈ 1 − np.

Ici, Y = nX̄n suit la loi binomiale de paramètres (n, p) et l’approximation normale est bonne
seulement si n et np sont grands.

Dr
 56

aft
Autre point de vue: Si p est très petit (np est petit), on a de fortes chances d’avoir
X1 = · · · = Xn = 0, ce qui donne X̄n = Sn2 = 0 et un intervalle de confiance de largeur 0.

P[X̄n = Sn2 = 0] = (P[Xi = 0])n = (1 − p)n ≈ 1 − np.

Ici, Y = nX̄n suit la loi binomiale de paramètres (n, p) et l’approximation normale est bonne
seulement si n et np sont grands.

Dr
Si n est grand et np est petit, alors Y approx. Poisson(np).
 56

aft
Autre point de vue: Si p est très petit (np est petit), on a de fortes chances d’avoir
X1 = · · · = Xn = 0, ce qui donne X̄n = Sn2 = 0 et un intervalle de confiance de largeur 0.

P[X̄n = Sn2 = 0] = (P[Xi = 0])n = (1 − p)n ≈ 1 − np.

Ici, Y = nX̄n suit la loi binomiale de paramètres (n, p) et l’approximation normale est bonne
seulement si n et np sont grands.

Dr
Si n est grand et np est petit, alors Y approx. Poisson(np).

Il existe des techniques de simulation spéciales pour ce contexte d’événements rares.

 57

Choix de la loi d’échantillonnage: “importance sampling”

Supposons que Y a la densité π sur R et qu’on veut estimer

aft
Z ∞
µ = Eπ [h(Y )] = h(y )π(y )dy < ∞,
−∞

où h : R → [0, ∞). Monte Carlo: générer Y ∼ π, retourner h(Y ).

Soit g une autre densité sur R, avec g (y ) > 0 quand h(y )π(y ) > 0. On a
Z ∞   
h(y )π(y ) h(Y )π(Y )
µ= g (y )dy = Eg = Eg [h(Y )L]
−∞ g (y ) g (Y )

Dr
où L = L(Y ) = π(Y )/g (Y ) est un rapport de vraisemblance.
 57

Choix de la loi d’échantillonnage: “importance sampling”

Supposons que Y a la densité π sur R et qu’on veut estimer

aft
Z ∞
µ = Eπ [h(Y )] = h(y )π(y )dy < ∞,
−∞

où h : R → [0, ∞). Monte Carlo: générer Y ∼ π, retourner h(Y ).

Soit g une autre densité sur R, avec g (y ) > 0 quand h(y )π(y ) > 0. On a
Z ∞   
h(y )π(y ) h(Y )π(Y )
µ= g (y )dy = Eg = Eg [h(Y )L]
−∞ g (y ) g (Y )

Dr
où L = L(Y ) = π(Y )/g (Y ) est un rapport de vraisemblance. Pour estimer µ (sans biais),
on peut générer Y1 , . . . , Yn i.i.d. selon g , et prendre la moyenne:

µ̂n =
1X
n
n
h(Yi )
i=1
π(Yi )
g (Yi )
=
1X
n
n
h(Yi )Li .
i=1

Par exemple, g pourrait être une densité normale même si π ne l’est pas.
 58

On a Var[µ̂n ] = Varg [h(Y )L]/n, où

= Eg h2 (Y )L2 − µ2
 
Varg [h(Y )L]

aft
Z ∞ 2 Z ∞ 2
h (y )π 2 (y ) 2 h (y )π 2 (y )
= 2
g (y )dy − µ = dy − µ2 .
−∞ g (y ) −∞ g (y )

Si g (y ) ≪ h(y )π(y ) quelque part, ou si h2 (y )π(y )/g (y ) converge vers 0 trop lentement
lorsque y → ±∞, alors la variance peut être énorme, voire infinie. Le choix de g est crucial.

Dr
 58

On a Var[µ̂n ] = Varg [h(Y )L]/n, où

= Eg h2 (Y )L2 − µ2
 
Varg [h(Y )L]

aft
Z ∞ 2 Z ∞ 2
h (y )π 2 (y ) 2 h (y )π 2 (y )
= 2
g (y )dy − µ = dy − µ2 .
−∞ g (y ) −∞ g (y )

Si g (y ) ≪ h(y )π(y ) quelque part, ou si h2 (y )π(y )/g (y ) converge vers 0 trop lentement
lorsque y → ±∞, alors la variance peut être énorme, voire infinie. Le choix de g est crucial.
Si L ≤ 1 quand h(Y ) > 0, alors

Varg [h(Y )L] = Eg [h(Y )2 L2 ] − µ2 = Eπ [h(Y )2 L] − µ2 ≤ Eπ [h(Y )2 ] − µ2 = Var[X ].

Dr
 58

On a Var[µ̂n ] = Varg [h(Y )L]/n, où

= Eg h2 (Y )L2 − µ2
 
Varg [h(Y )L]

aft
Z ∞ 2 Z ∞ 2
h (y )π 2 (y ) 2 h (y )π 2 (y )
= 2
g (y )dy − µ = dy − µ2 .
−∞ g (y ) −∞ g (y )

Si g (y ) ≪ h(y )π(y ) quelque part, ou si h2 (y )π(y )/g (y ) converge vers 0 trop lentement
lorsque y → ±∞, alors la variance peut être énorme, voire infinie. Le choix de g est crucial.
Si L ≤ 1 quand h(Y ) > 0, alors

Varg [h(Y )L] = Eg [h(Y )2 L2 ] − µ2 = Eπ [h(Y )2 L] − µ2 ≤ Eπ [h(Y )2 ] − µ2 = Var[X ].

Dr
Quel est le g optimal? Si on prend g (y ) proportionnel à h(y )π(y ), alors l’estimateur IS
devient une constante, il ne dépend plus de y , donc sa variance estRzéro! Comme
∞
l’estimateur est sans biais, cette constante est nécessairement µ = −∞ h(y )π(y )dy .
Implanter ceci exactement est rarement possible, mais on peut souvent l’approximer.
 59

Mini-Exemple: estimer une probabilité

On veut estimer p = P[Y > y0 ] = E[I[Y > y0 ]], où Y est exponentielle de paramètre (taux)

aft
λ, i.e., π(y ) = λe −λy pour y ≥ 0.
(Exemple purement académique; on sait que P[Y > y ] = e −λy .)

λ .....
.....
.....
.....
.....
.....
......
......
......
......
.......
.......
.......
λ/2 ........
........
.........
.........
..........
...........
............
.............. ↙ f (y )
................
.......................

0
0 1/λ

Monte Carlo: générer Y selon la bonne densité exponentielle π, soit Y = − ln(1 − U)/λ, et
calculer X = I[Y > y0 ]. Répéter n fois et calculer la moyenne des n réalisations de X .
Dr ..........................

2/λ y0
...........................................
..........................................

y0
 60

Mini-Exemple: estimer une probabilité

On veut estimer p = P[Y > y0 ] = E[I[Y > y0 ]], où Y est exponentielle de paramètre (taux)

aft
λ, i.e., π(y ) = λe −λy pour y ≥ 0.

Soit g = π0 une autre densité exponentielle, de paramètre λ0 ̸= λ. On a

Z ∞ Z ∞
π(y )
p= I[y ≥ y0 ]π(y )dy = I[y ≥ y0 ] π0 (y )dy = Eπ0 [Xis ],
0 0 π0 (y )

où

Xis = I[Y0 ≥ y0 ]

Dr
π(Y0 )
π0 (Y0 )
= I[Y0 ≥ y0 ]
λ exp[−λY0 ]
λ0 exp[−λ0 Y0 ]

et Y0 = − ln(1 − U)/λ0 est une v.a. exponentielle de taux λ0 .

Cela donne un estimateur sans biais peu importe λ0 > 0.
λ
= I[Y0 ≥ y0 ] exp[−(λ − λ0 )Y0 ]
λ0
 61
On a aussi

Var[Xis ] = E[Xis2 ] − p 2

aft
Z ∞ 2
π (y )
= 2
π0 (y )dy − p 2
y0 π0 (y )
Z ∞
= (λ/λ0 )2 exp[−2(λ − λ0 )y ]λ0 exp[−λ0 y ]dy − p 2
y
0
λ2
exp[−(2λ − λ0 )y0 ] − p 2 si 0 < λ0 < 2λ,

= λ0 (2λ − λ0 )
∞ sinon.


Dr
Pour λ0 = λ, la variance est p(1 − p). Lorsque λ0 s’approche de 0 ou de 2λ, la variance tend
vers l’infini.
 61
On a aussi

Var[Xis ] = E[Xis2 ] − p 2

aft
Z ∞ 2
π (y )
= 2
π0 (y )dy − p 2
y0 π0 (y )
Z ∞
= (λ/λ0 )2 exp[−2(λ − λ0 )y ]λ0 exp[−λ0 y ]dy − p 2
y
0
λ2
exp[−(2λ − λ0 )y0 ] − p 2 si 0 < λ0 < 2λ,

= λ0 (2λ − λ0 )
∞ sinon.


Dr
Pour λ0 = λ, la variance est p(1 − p). Lorsque λ0 s’approche de 0 ou de 2λ, la variance tend
vers l’infini. Exercice: La variance est minimisée pour λ0 = λ + 1/y0 − (λ2 + 1/y02 )1/2 < λ.
 61
On a aussi

Var[Xis ] = E[Xis2 ] − p 2

aft
Z ∞ 2
π (y )
= 2
π0 (y )dy − p 2
y0 π0 (y )
Z ∞
= (λ/λ0 )2 exp[−2(λ − λ0 )y ]λ0 exp[−λ0 y ]dy − p 2
y
0
λ2
exp[−(2λ − λ0 )y0 ] − p 2 si 0 < λ0 < 2λ,

= λ0 (2λ − λ0 )
∞ sinon.


Dr
Pour λ0 = λ, la variance est p(1 − p). Lorsque λ0 s’approche de 0 ou de 2λ, la variance tend
vers l’infini. Exercice: La variance est minimisée pour λ0 = λ + 1/y0 − (λ2 + 1/y02 )1/2 < λ.

Cet exemple montre que remplacer π par une autre densité g peut réduire la variance, mais
peut aussi l’augmenter, et même la rendre infinie.
La variance est souvent très sensible au choix de g .
 62
Exemple numérique: λ = 1 et y0 = 4.

Var[Xis ]/Var[X ]

aft
2

1.5

0.5

Dr
0 0.25 0.5 0.75 1

Variance min. pour λ0 ≈ 0.2192; donne Var[Xis ]/Var[X ] ≈ 0.0962.

1.25
λ0

Ainsi, pour estimer p = P[Y > 4], on change la moyenne de Y pour 1/λ0 ≈ 1/.2192 ≈ 4.56,
qui est proche de 4. Semble raisonnable.
 63
Densité g optimale pour cet exemple: g (y ) ∝ I[y > y0 ]π(y ).
C’est une exponentielle tronquée à [y0 , ∞):
π(y ) λe −λy

aft
g (y ) = = −λy = λe −λ(y −y0 ) pour y > y0 , et 0 ailleurs.
P[Y > y0 ] e 0
C’est la densité de y0 + X où X est une v.a. exponentielle de taux λ.
On peut donc générer Y = y0 − ln(1 − U)/λ.

Dr
 63
Densité g optimale pour cet exemple: g (y ) ∝ I[y > y0 ]π(y ).
C’est une exponentielle tronquée à [y0 , ∞):
π(y ) λe −λy

aft
g (y ) = = −λy = λe −λ(y −y0 ) pour y > y0 , et 0 ailleurs.
P[Y > y0 ] e 0
C’est la densité de y0 + X où X est une v.a. exponentielle de taux λ.
On peut donc générer Y = y0 − ln(1 − U)/λ.
Le rapport de vraisemblance est L = π(Y )/g (Y ) = e −λy0 et on a

Xis = I[Y ≥ y0 ] exp[−λy0 ] = exp[−λy0 ] = p,

Dr
car on a toujours Y ≥ y0 .
Ce Xis est donc un estimateur sans biais de variance zero!
 63
Densité g optimale pour cet exemple: g (y ) ∝ I[y > y0 ]π(y ).
C’est une exponentielle tronquée à [y0 , ∞):
π(y ) λe −λy

aft
g (y ) = = −λy = λe −λ(y −y0 ) pour y > y0 , et 0 ailleurs.
P[Y > y0 ] e 0
C’est la densité de y0 + X où X est une v.a. exponentielle de taux λ.
On peut donc générer Y = y0 − ln(1 − U)/λ.
Le rapport de vraisemblance est L = π(Y )/g (Y ) = e −λy0 et on a

Xis = I[Y ≥ y0 ] exp[−λy0 ] = exp[−λy0 ] = p,

car on a toujours Y ≥ y0 .

Dr
Ce Xis est donc un estimateur sans biais de variance zero!

En théorie, de tels estimateurs “magiques” existent toujours. En pratique, il sont très

difficiles à trouver et implanter. Par contre, on peut souvent les approximer et construire
ainsi des estimateurs plus performants.
Rendu ici, 8 novembre 2024
 64
Exemple. Soient Y1 et Y2 des v.a. indép. de densités π1 et π2 , sur R.
(
Y1 + Y2 − K si Y1 ≤ a et Y1 + Y2 ≥ b,
X =

aft
0 sinon,

où K > 0, et a et b sont des constantes. On veut estimer µ = E[X ].

Dr
 64
Exemple. Soient Y1 et Y2 des v.a. indép. de densités π1 et π2 , sur R.
(
Y1 + Y2 − K si Y1 ≤ a et Y1 + Y2 ≥ b,
X =

aft
0 sinon,

où K > 0, et a et b sont des constantes. On veut estimer µ = E[X ].

MC standard: générer Y1 et Y2 selon π1 et π2 , et calculer X .

Dr
 64
Exemple. Soient Y1 et Y2 des v.a. indép. de densités π1 et π2 , sur R.
(
Y1 + Y2 − K si Y1 ≤ a et Y1 + Y2 ≥ b,
X =

aft
0 sinon,

où K > 0, et a et b sont des constantes. On veut estimer µ = E[X ].

MC standard: générer Y1 et Y2 selon π1 et π2 , et calculer X .

Stratégie IS: Éviter de gaspiller des échantillons dans la région où X = 0.
Par exemple, générer Y1 selon sa densité conditionnelle à Y1 ≤ a,
puis générer Y2 selon sa densité conditionnelle à Y1 + Y2 ≥ b,

Dr
i.e., tronquée à l’intervalle [b − Y1 , ∞).
 64
Exemple. Soient Y1 et Y2 des v.a. indép. de densités π1 et π2 , sur R.
(
Y1 + Y2 − K si Y1 ≤ a et Y1 + Y2 ≥ b,
X =

aft
0 sinon,

où K > 0, et a et b sont des constantes. On veut estimer µ = E[X ].

MC standard: générer Y1 et Y2 selon π1 et π2 , et calculer X .

Stratégie IS: Éviter de gaspiller des échantillons dans la région où X = 0.
Par exemple, générer Y1 selon sa densité conditionnelle à Y1 ≤ a,
puis générer Y2 selon sa densité conditionnelle à Y1 + Y2 ≥ b,

Dr
i.e., tronquée à l’intervalle [b − Y1 , ∞).
La nouvelle densité de Y1 est

g1 (y ) =
π1 (y )
P[Y1 ≤ a]

et celle de Y2 conditionelle à Y1 = y1 est

=
π1 (y )
F1 (a)
pour y ≤ a,
 65

π2 (y ) π2 (y )
g2 (y | y1 ) = = pour y ≥ b − y1 ,

aft
P[Y2 ≥ b − y1 ] 1 − F2 (b − y1 )
où F1 et F2 sont les fonctions de répartition de Y1 et Y2 .

Dr
 65

π2 (y ) π2 (y )
g2 (y | y1 ) = = pour y ≥ b − y1 ,

aft
P[Y2 ≥ b − y1 ] 1 − F2 (b − y1 )
où F1 et F2 sont les fonctions de répartition de Y1 et Y2 . On a
Z ∞Z ∞
µ = X π2 (y2 )π1 (y1 )dy2 dy1
−∞ −∞
Z a Z ∞
π2 (y2 )π1 (y1 )
= X g2 (y2 | y1 )g1 (y1 )dy2 dy1
−∞ b−y1 g2 (y2 | y1 )g1 (y1 )
Z a Z ∞
=

Dr
−∞
= Eg [Xis ]
X F1 (a) (1 − F2 (b − y1 ))g2 (y2 | y1 )g1 (y1 )dy2 dy1
b−y1

où Xis = X L avec L = F1 (a) (1 − F2 (b − Y1 )), et Eg désigne l’espérance sous g1 et g2 .

On a toujours Var[Xis ] = Eg [X 2 L2 ] − µ2 = Eπ [X 2 L] − µ2 < Var[X ] car L < 1.
 66

Expérience: Supposons π1 et π2 ∼ N (1, 1), K = 1, b = 2, et a = 1/2.

Ici, Y1 − 1 et Y2 − 1 sont N (0, 1), de cdf Φ.

aft
Essayer MC vs IS, avec n = 105 , et comparer les variances.

Dr
 66

Expérience: Supposons π1 et π2 ∼ N (1, 1), K = 1, b = 2, et a = 1/2.

Ici, Y1 − 1 et Y2 − 1 sont N (0, 1), de cdf Φ.

aft
Essayer MC vs IS, avec n = 105 , et comparer les variances.

On a F1 (a) = P[Y1 < a] = P[Y1 − 1 < a − 1] = Φ(a − 1).

On pose U1 ∼ U(0, Φ(a − 1)) et Y1 = 1 + Φ−1 (U1 ).

0Dr
Φ(a − 1)
U1

Y1 − 1 a−1
Φ(y )

y
 67
Ensuite, pour Y1 = y1 donné, on veut Y2 > b − y1 . On a
1 − F2 (b − y1 ) = P[Y2 > b − y1 ] = P[Y2 − 1 > b − 1 − y1 ] = 1 − Φ(b − 1 − y1 ).
On pose U2 ∼ U(Φ(b − 1 − y1 ), 1)) et Y2 = 1 + Φ−1 (U2 ).

aft
U2 1
Φ(b1 − 1 − y1 )

Φ(y )

0 y
0 b1 − 1 − y1 Y2 − 1

Dr
 67
Ensuite, pour Y1 = y1 donné, on veut Y2 > b − y1 . On a
1 − F2 (b − y1 ) = P[Y2 > b − y1 ] = P[Y2 − 1 > b − 1 − y1 ] = 1 − Φ(b − 1 − y1 ).
On pose U2 ∼ U(Φ(b − 1 − y1 ), 1)) et Y2 = 1 + Φ−1 (U2 ).

aft
U2 1
Φ(b1 − 1 − y1 )

Φ(y )

0 y
0 b1 − 1 − y1 Y2 − 1

Dr
On calcule l’estimateur Xis = X Φ(a − 1)(1 − Φ(b − 1 − Y1 )).
 67
Ensuite, pour Y1 = y1 donné, on veut Y2 > b − y1 . On a
1 − F2 (b − y1 ) = P[Y2 > b − y1 ] = P[Y2 − 1 > b − 1 − y1 ] = 1 − Φ(b − 1 − y1 ).
On pose U2 ∼ U(Φ(b − 1 − y1 ), 1)) et Y2 = 1 + Φ−1 (U2 ).

aft
U2 1
Φ(b1 − 1 − y1 )

Φ(y )

0 y
0 b1 − 1 − y1 Y2 − 1

Dr
On calcule l’estimateur Xis = X Φ(a − 1)(1 − Φ(b − 1 − Y1 )).
La variance Sn2 est environ 40 fois plus petite avec Xis qu’avec X .

Estimator
X
Xis
µ̂n
0.0733
0.0742
Sn2
0.1188
0.0027
IC à 95%
(0.071, 0.075)
(0.074, 0.075)
 68

Exemple: Une chaine de Markov

{Xn , n ≥ 0} sur les états {0, 1, . . . , K }, avec X0 = x0 .
Probabilités de transition Pi,j = P[Xn = j | Xn−1 = i].

aft
On veut estimer la probabilité µ = µ(x0 ) d’atteindre K avant de revenir à 0.

Dr
 68

Exemple: Une chaine de Markov

{Xn , n ≥ 0} sur les états {0, 1, . . . , K }, avec X0 = x0 .
Probabilités de transition Pi,j = P[Xn = j | Xn−1 = i].

aft
On veut estimer la probabilité µ = µ(x0 ) d’atteindre K avant de revenir à 0.
Estimateur naif: I[XT = K ] où T = inf{n ≥ 1 : Xn ∈ {0, K }}.

Dr
 68

Exemple: Une chaine de Markov

{Xn , n ≥ 0} sur les états {0, 1, . . . , K }, avec X0 = x0 .
Probabilités de transition Pi,j = P[Xn = j | Xn−1 = i].

aft
On veut estimer la probabilité µ = µ(x0 ) d’atteindre K avant de revenir à 0.
Estimateur naif: I[XT = K ] où T = inf{n ≥ 1 : Xn ∈ {0, K }}.
IS: changer les Pi,j pour des Qi,j pour augmenter les chances d’aller à K .

Dr
 68

Exemple: Une chaine de Markov

{Xn , n ≥ 0} sur les états {0, 1, . . . , K }, avec X0 = x0 .
Probabilités de transition Pi,j = P[Xn = j | Xn−1 = i].

aft
On veut estimer la probabilité µ = µ(x0 ) d’atteindre K avant de revenir à 0.
Estimateur naif: I[XT = K ] où T = inf{n ≥ 1 : Xn ∈ {0, K }}.
IS: changer les Pi,j pour des Qi,j pour augmenter les chances d’aller à K .
QT
La probabilité d’une trajectoire X1 , X2 , . . . , XT est n=1 PXn−1 ,Xn avec les probabilités originales et
QT
n=1 QXn−1 ,Xn avec les nouvelles probabilités. En faisant la somme sur toutes les valeurs de T et
toutes les trajectoires possibles pour chaque T , on obtient
T T T
µ =

=
X

T ,X1 ,...,XT

T ,X1 ,...,XT
Dr
I[XT = K ]

I[XT = K ]L
Y

n=1
T
Y

n=1
PXn−1 ,Xn =
X

T ,X1 ,...,XT
I[XT = K ]

QXn−1 ,Xn = EQ [I[XT = K ]L] = EQ [Xis ].

De plus, si on a toujours L ≤ 1 lorsque XT = K , alors VarQ [I[XT = K ]L] ≤ VarP [I[XT = K ].

Y PXn−1 ,Xn Y
Q
n=1 Xn−1 ,Xn n=1
QXn−1 ,Xn
 69

aft
Comment choisir les Qi,j ?
Idée simpliste: bloquer les retours à 0 en posant Qi,0 = 0 pour tout i > 0, et renormaliser les
autres probabilités: Qi,j = Pi,j /(1 − Pi,0 ) pour i, j > 0, et Q0,j = P0,j pour tout j.

Dr
 69

aft
Comment choisir les Qi,j ?
Idée simpliste: bloquer les retours à 0 en posant Qi,0 = 0 pour tout i > 0, et renormaliser les
autres probabilités: Qi,j = Pi,j /(1 − Pi,0 ) pour i, j > 0, et Q0,j = P0,j pour tout j.
On a alors P[XT = K ] = 1 et le nouvel estimateur est
T T
Y PXn−1 ,Xn Y
L = L(X1 , . . . , XT ) = = (1 − PXn−1 ,0 ) < 1.
QXn−1 ,Xn
n=1 n=2

La variance est toujours réduite car L < 1 quand XT = K .

Dr
Mais: chaque simulation risque d’être très longue (on peut de se promener autour de l’état 0
très longtemps). L’efficacité n’est pas nécessairement améliorée. On doit faire mieux.
 70
Une marche aléatoire unidimensionnelle sur {0, 1, . . . , K }, avec x0 = 1.
Soient Pi,i+1 = p et Pi,i−1 = 1 − p pour 1 ≤ i ≤ K − 1, et P0,1 = PK ,K −1 = 1.
p p

aft
1 p p p

0 1 2 3 ··· K −1 K
1−p 1−p 1−p 1−p 1−p 1

Si p < 1/2 et K est grand, la chaine est attirée vers 0.

Couper l’accès à 0 ne suffit pas, il faut augmenter l’attirance vers K .

Dr
 70
Une marche aléatoire unidimensionnelle sur {0, 1, . . . , K }, avec x0 = 1.
Soient Pi,i+1 = p et Pi,i−1 = 1 − p pour 1 ≤ i ≤ K − 1, et P0,1 = PK ,K −1 = 1.
p p

aft
1 p p p

0 1 2 3 ··· K −1 K
1−p 1−p 1−p 1−p 1−p 1

Si p < 1/2 et K est grand, la chaine est attirée vers 0.

Couper l’accès à 0 ne suffit pas, il faut augmenter l’attirance vers K .
On va modifier les probabilités Pi,j pour Qi,i+1 = q et Qi,i−1 = 1 − q, pour q > p.

Dr
 70
Une marche aléatoire unidimensionnelle sur {0, 1, . . . , K }, avec x0 = 1.
Soient Pi,i+1 = p et Pi,i−1 = 1 − p pour 1 ≤ i ≤ K − 1, et P0,1 = PK ,K −1 = 1.
p p

aft
1 p p p

0 1 2 3 ··· K −1 K
1−p 1−p 1−p 1−p 1−p 1

Si p < 1/2 et K est grand, la chaine est attirée vers 0.

Couper l’accès à 0 ne suffit pas, il faut augmenter l’attirance vers K .
On va modifier les probabilités Pi,j pour Qi,i+1 = q et Qi,i−1 = 1 − q, pour q > p.

L =
Y
Dr
Examinons le rapport de vraisemblance lorsque XT = K . Pour chaque paire d’états (i, i + 1)
on va de i à i + 1 une fois de plus que de i + 1 à i. Donc
T
PXn−1 ,Xn
QXn−1 ,Xn
n=1
=
 K −1 
p
q
p(1 − p) (T −K )/2
q(1 − q)

.
 71

 K −1 
p(1 − p) (T −K )/2

p
L = .
q(1 − q)

aft
q
Pour s’assurer que L < 1, prenons q > p et q(1 − q) ≥ p(1 − p), i.e., p < q ≤ 1 − p.

Dr
 71

 K −1 
p(1 − p) (T −K )/2

p
L = .
q(1 − q)

aft
q
Pour s’assurer que L < 1, prenons q > p et q(1 − q) ≥ p(1 − p), i.e., p < q ≤ 1 − p.
En maximisant q sous cette contrainte (pour aller à K le plus vite et le plus souvent
possible), on obtient q = 1 − p. Le RV se simplifie alors et devient une constante:
 K −1  K −1
p p
L= = < 1.
q 1−p

Dr
Avant, l’estimateur était 1 avec prob. µ, et 0 sinon. Variance de µ(1 − µ).
Maintenant il est L avec prob. µ/L, et 0 sinon.
Nouvelle variance: L2 µ/L − µ2 = µ(L − µ) < Lµ(1 − µ).
Avec IS, la variance est inférieure à L fois l’ancienne variance.
 71

 K −1 
p(1 − p) (T −K )/2

p
L = .
q(1 − q)

aft
q
Pour s’assurer que L < 1, prenons q > p et q(1 − q) ≥ p(1 − p), i.e., p < q ≤ 1 − p.
En maximisant q sous cette contrainte (pour aller à K le plus vite et le plus souvent
possible), on obtient q = 1 − p. Le RV se simplifie alors et devient une constante:
 K −1  K −1
p p
L= = < 1.
q 1−p

Dr
Avant, l’estimateur était 1 avec prob. µ, et 0 sinon. Variance de µ(1 − µ).
Maintenant il est L avec prob. µ/L, et 0 sinon.
Nouvelle variance: L2 µ/L − µ2 = µ(L − µ) < Lµ(1 − µ).
Avec IS, la variance est inférieure à L fois l’ancienne variance.

Example: si p = 1/3 et K = 101, on a L = 2−100 ≈ 0.833333 × 10−30 .

Dans ce cas, la variance est divisée par presque 10−30 !
 72

Monte Carlo Conditionnel (conditional Monte Carlo (CMC))

aft
D’abord un exemple très simple.
Supposons que l’on veut estimer µ = P[Y1 + · · · + Yt > x] où Y1 , . . . , Yt sont indépendants.
Estimateur MC évident: X = I[Y1 + · · · + Yt > x].

Dr
 72

Monte Carlo Conditionnel (conditional Monte Carlo (CMC))

aft
D’abord un exemple très simple.
Supposons que l’on veut estimer µ = P[Y1 + · · · + Yt > x] où Y1 , . . . , Yt sont indépendants.
Estimateur MC évident: X = I[Y1 + · · · + Yt > x].
Estimateur CMC: on ne génère pas Yt . On estime µ par Xe,t−1 = P[Y1 + · · · + Yt > x |
Y1 , . . . , Yt−1 ] = P[Yt > x − Y1 − · · · − Yt−1 | Y1 , . . . , Yt−1 ] = 1 − Ft [x − Y1 − · · · − Yt−1 ].

Dr
 72

Monte Carlo Conditionnel (conditional Monte Carlo (CMC))

aft
D’abord un exemple très simple.
Supposons que l’on veut estimer µ = P[Y1 + · · · + Yt > x] où Y1 , . . . , Yt sont indépendants.
Estimateur MC évident: X = I[Y1 + · · · + Yt > x].
Estimateur CMC: on ne génère pas Yt . On estime µ par Xe,t−1 = P[Y1 + · · · + Yt > x |
Y1 , . . . , Yt−1 ] = P[Yt > x − Y1 − · · · − Yt−1 | Y1 , . . . , Yt−1 ] = 1 − Ft [x − Y1 − · · · − Yt−1 ].
Estimateur CMC plus général: Xe,s = P[Y1 + · · · + Yt > x | Y1 , . . . , Ys ], pour s ≤ t.

Dr
 72

Monte Carlo Conditionnel (conditional Monte Carlo (CMC))

aft
D’abord un exemple très simple.
Supposons que l’on veut estimer µ = P[Y1 + · · · + Yt > x] où Y1 , . . . , Yt sont indépendants.
Estimateur MC évident: X = I[Y1 + · · · + Yt > x].
Estimateur CMC: on ne génère pas Yt . On estime µ par Xe,t−1 = P[Y1 + · · · + Yt > x |
Y1 , . . . , Yt−1 ] = P[Yt > x − Y1 − · · · − Yt−1 | Y1 , . . . , Yt−1 ] = 1 − Ft [x − Y1 − · · · − Yt−1 ].
Estimateur CMC plus général: Xe,s = P[Y1 + · · · + Yt > x | Y1 , . . . , Ys ], pour s ≤ t.
Pour s = t: aucun changement.

Dr
Pour s = t − 1, c’est le cas précédent.
Pour s = 0: Xe,0 = P[Y1 + · · · + Yt > x] = µ (variance réduite à zéro).
Plus s est petit, plus la variance est réduite, mais plus l’estimateur est difficile à calculer.
 73

Cadre général
Idée: remplacer l’estimateur X par E[X | Z ] où Z est une autre v.a., ou plus généralement

aft
par E[X | G], où G représente une information partielle sur X .

Dr
 73

Cadre général
Idée: remplacer l’estimateur X par E[X | Z ] où Z est une autre v.a., ou plus généralement

aft
par E[X | G], où G représente une information partielle sur X . L’estimateur CMC est
def
Xe = E[X | G].

On a
E[Xe ] = E[E[X | G]] = E[X ].

Dr
 73

Cadre général
Idée: remplacer l’estimateur X par E[X | Z ] où Z est une autre v.a., ou plus généralement

aft
par E[X | G], où G représente une information partielle sur X . L’estimateur CMC est
def
Xe = E[X | G].

On a
E[Xe ] = E[E[X | G]] = E[X ].
De plus,

Var[X ] = E[ Var[X | G] ] + Var[E[X | G]] = E[Var[X | G]] + Var[Xe ],

et donc
Dr
| {z }

pour G connu
|
Var. résiduelle

(éliminé par CMC)

{z }
Var. due à la
variation de G

Var[Xe ] = Var[X ] − E[Var[X | G]] ≤ Var[X ].

 74

Pour minimiser la variance, on doit maximiser E[Var[X | G]], i.e., G doit contenir le moins
d’information possible.

aft
On sait en effet que G1 ⊂ G2 implique que E[Var[X | G1 ]] ≥ E[Var[X | G2 ]]. Preuve:

Var[X | G1 ] = E[Var[X | G2 ] | G1 ] + Var[E[X | G2 ] | G1 ],

et il suffit de prendre l’espérance des deux cotés:

E[Var[X | G1 ]] ≥ E[E[Var[X | G2 ] | G1 ]] = E[Var[X | G2 ]].

Dr
 74

Pour minimiser la variance, on doit maximiser E[Var[X | G]], i.e., G doit contenir le moins
d’information possible.

aft
On sait en effet que G1 ⊂ G2 implique que E[Var[X | G1 ]] ≥ E[Var[X | G2 ]]. Preuve:

Var[X | G1 ] = E[Var[X | G2 ] | G1 ] + Var[E[X | G2 ] | G1 ],

et il suffit de prendre l’espérance des deux cotés:

E[Var[X | G1 ]] ≥ E[E[Var[X | G2 ] | G1 ]] = E[Var[X | G2 ]].

Dr
Mais moins G1 contient d’information, plus il est difficile de calculer Xe .
On doit donc faire un compromis.
Dans certains cas, Xe peut être moins coûteux à calculer que X .
Cas limites:
Si G ne contient aucune information pertinente à X : Xe = E[X | G] = E[X ] = µ.
Si G permet de calculer X (i.e., X est G-mesurable): Xe = X .
 75

Déplacement d’une poutre en porte-à-faux

Déplacement X pour une charge horizontale Y2 et charge verticale Y3 :

aft
s
κ Y2 2 Y3 2
X = h(Y1 , Y2 , Y3 ) = + 4
Y1 w 4 t
où Y1 , Y2 , Y3 sont normales indépendantes, Yj ∼ N (µj , σj2 ).
Supposons que l’on veut estimer µ = P[X ≤ x] pour x fixé.
(a) MC ordinaire: générer X et calculer I = I[X ≤ x], qui vaut 0 ou 1.
(b) CMC: générer seulement Y2 , Y3 , et calculer J = P[X ≤ x | Y2 , Y3 ]. On a

Donc pour x > 0, Dr

X ≤ x ssi Y1 ≥
r
κ Y22 Y32 def
x w4
+ 4 = W1 (x).
t

J = P[X ≤ x | Y2 , Y3 ] = P[Y1 ≥ W1 (x) | W1 (x)] = 1 − Φ((W1 (x) − µ1 )/σ1 ).

 75

Déplacement d’une poutre en porte-à-faux

Déplacement X pour une charge horizontale Y2 et charge verticale Y3 :

aft
s
κ Y2 2 Y3 2
X = h(Y1 , Y2 , Y3 ) = + 4
Y1 w 4 t
où Y1 , Y2 , Y3 sont normales indépendantes, Yj ∼ N (µj , σj2 ).
Supposons que l’on veut estimer µ = P[X ≤ x] pour x fixé.
(a) MC ordinaire: générer X et calculer I = I[X ≤ x], qui vaut 0 ou 1.
(b) CMC: générer seulement Y2 , Y3 , et calculer J = P[X ≤ x | Y2 , Y3 ]. On a

Donc pour x > 0, Dr

X ≤ x ssi Y1 ≥
r
κ Y22 Y32 def
x w4
+ 4 = W1 (x).
t

J = P[X ≤ x | Y2 , Y3 ] = P[Y1 ≥ W1 (x) | W1 (x)] = 1 − Φ((W1 (x) − µ1 )/σ1 ).

Dans les deux cas, pour I et pour J, on génère n réalisations indépendantes de l’estimateur
puis on calcule la moyenne et un intervalle de confiance. Suite dans le devoir 7.
 76

Exemple: Réseau d’activités stochastique

On veut estimer µ = P[T > x]. Estimateur naı̈f: X = I[T > x], vaut 0 ou 1.

aft
Dr
 76

Exemple: Réseau d’activités stochastique

On veut estimer µ = P[T > x]. Estimateur naı̈f: X = I[T > x], vaut 0 ou 1.

aft
Soit L ⊆ A un ensemble d’activités (arcs) tel que chaque chemin de la source au puits
contient exactement un arc de L. (L est une coupe orientée.)
Soit G = {Yj , j ̸∈ L}. Exemple: L = {4, 5, 6, 8, 9}.

5 8 sink
Y10
Y5
Y9 Y12

Dr
source 0
Y1

Y0
Y2
2

1
Y3
Y4
4

3
Y7
Y8

Y6
7

6
Y11
 77

Estimateur CMC:
Xe = P[T > x | G] = P[T > x | {Yj , j ̸∈ L}] = P[T > x | Y0 , Y1 , Y2 , Y3 , Y7 , Y10 , Y11 , Y12 ].

aft
On le calcule comme suit.
Pour chaque l ∈ L, allant disons de al à bl , calculer la longueur αl du plus long chemin de la
source à al , puis la longueur βl du plus long chemin de bl au puits.
Aucun chemin passant par l est plus long que x ssi αl + Yl + βl ≤ x.

Dr
 77

Estimateur CMC:
Xe = P[T > x | G] = P[T > x | {Yj , j ̸∈ L}] = P[T > x | Y0 , Y1 , Y2 , Y3 , Y7 , Y10 , Y11 , Y12 ].

aft
On le calcule comme suit.
Pour chaque l ∈ L, allant disons de al à bl , calculer la longueur αl du plus long chemin de la
source à al , puis la longueur βl du plus long chemin de bl au puits.
Aucun chemin passant par l est plus long que x ssi αl + Yl + βl ≤ x.
Conditionnellement à G, on a cette condition avec probabilité
P[Yl ≤ x − αl − βl ] = Fl [x − αl − βl ].

Dr
 77

Estimateur CMC:
Xe = P[T > x | G] = P[T > x | {Yj , j ̸∈ L}] = P[T > x | Y0 , Y1 , Y2 , Y3 , Y7 , Y10 , Y11 , Y12 ].

aft
On le calcule comme suit.
Pour chaque l ∈ L, allant disons de al à bl , calculer la longueur αl du plus long chemin de la
source à al , puis la longueur βl du plus long chemin de bl au puits.
Aucun chemin passant par l est plus long que x ssi αl + Yl + βl ≤ x.
Conditionnellement à G, on a cette condition avec probabilité
P[Yl ≤ x − αl − βl ] = Fl [x − αl − βl ].

Dr
Puisque les Yl sont indépendants, on obtient

Xe = 1 − P[Yl ≤ x − αl − βl pour tout l] = 1 −

Y

l∈L
Fl [x − αl − βl ].

Pour l’exemple numérique donné en page 5, cela divise la variance environ par 4.
Et cet estimateur peut être moins coûteux à calculer que X , car moins de Yj ’s à générer.
Rendu ici, 11 nov. 2024
 78

Valeurs aléatoires communes (common random numbers (CRN))

Idée: pour comparer deux (ou plusieurs) systèmes semblables, utiliser les mêmes nombres

aft
aléatoires uniformes aux mêmes endroits pour tous les systèmes.
Supposons que µ1 = E[X1 ] et µ2 = E[X2 ]. On veut estimer µ2 − µ1 = E[X2 − X1 ].

Dr
 78

Valeurs aléatoires communes (common random numbers (CRN))

Idée: pour comparer deux (ou plusieurs) systèmes semblables, utiliser les mêmes nombres

aft
aléatoires uniformes aux mêmes endroits pour tous les systèmes.
Supposons que µ1 = E[X1 ] et µ2 = E[X2 ]. On veut estimer µ2 − µ1 = E[X2 − X1 ].

On simulant X1 et X2 avec les mêmes nombres aléatoires, on ne change pas leurs lois de
probabilité individuelles, mais on peut induire une covariance positive entre les deux. On a

Var[X2 − X1 ] = Var[X2 ] + Var[X1 ] − 2 Cov[X1 , X2 ].

Dr
 78

Valeurs aléatoires communes (common random numbers (CRN))

Idée: pour comparer deux (ou plusieurs) systèmes semblables, utiliser les mêmes nombres

aft
aléatoires uniformes aux mêmes endroits pour tous les systèmes.
Supposons que µ1 = E[X1 ] et µ2 = E[X2 ]. On veut estimer µ2 − µ1 = E[X2 − X1 ].

On simulant X1 et X2 avec les mêmes nombres aléatoires, on ne change pas leurs lois de
probabilité individuelles, mais on peut induire une covariance positive entre les deux. On a

Var[X2 − X1 ] = Var[X2 ] + Var[X1 ] − 2 Cov[X1 , X2 ].

Dr
“Independent random numbers”(IRN): X1 et X2 sont indépendants.
“Common random numbers” (CRN): On utilise l’inversion avec les mêmes uniformes Uj à la
même place pour simuler X1 et X2 . On espère avoir Cov[X1 , X2 ] > 0. Parfois, on peut le
prouver.
 79

Réseau d’activités stochastique

5 8 puit
Y10

aft
Exemple: Y5
Y9 Y12
On augmente E[Y2 ] de 7.0 à 10.0,
et E[Y4 ] de 16.5 à 18.5. 2 4 7
Y4 Y8
Y1
Impact sur la durée du projet? Y2 Y7 Y11

source 0 1 3 6
Y0 Y3 Y6
X1 = durée du projet selon les lois originales.

Dr
X2 = durée du projet avec les lois modifiées.
On veut étudier la loi de ∆ = X1 − X2 et estimer E[∆].
On suppose que Yj = Fj−1 (Uj ) pour X1 et Ỹj = F̃j−1 (Ũj ) pour X2 .
IRN: Les Ũj sont indépendants des Uj .
CRN: Ũj = Uj pour chaque j.
 80

Essayons n = 100, 000 répétitions pour chaque estimateur.

Avec IRN, les réalisations de ∆ vont de −223.22 à 280.92.

aft
moyenne = 1.326, variance = 967.
Intervalle de confiance à 95% pour E[∆]: (1.133, 1.519).
Avec CRN, ∆ va de 0 à 49.88.
moyenne = 1.528, variance = 9.1,
IC à 95% pour E[∆]: (1.510, 1.547).
CRN réduit la variance par un facteur d’environ 106.

Dr
 80

Essayons n = 100, 000 répétitions pour chaque estimateur.

Avec IRN, les réalisations de ∆ vont de −223.22 à 280.92.

aft
moyenne = 1.326, variance = 967.
Intervalle de confiance à 95% pour E[∆]: (1.133, 1.519).
Avec CRN, ∆ va de 0 à 49.88.
moyenne = 1.528, variance = 9.1,
IC à 95% pour E[∆]: (1.510, 1.547).
CRN réduit la variance par un facteur d’environ 106.

Dr
Avec CRN, 67,880 réalisations de ∆ sont 0 (parce que les Yj modifiés ne sont pas sur le plus
long chemin) et les 32,120 autres réalisations de ∆ sont toutes positives. Raison: quand on
augmente sa moyenne θj , Yj = −θj ln(1 − Uj ) (loi exponentielle) ne peut pas diminuer pour
Uj fixé, car − ln(1 − Uj ) > 0. Ainsi la longueur du plus long chemin ne peut pas diminuer.
Avec IRN, on a deux valeurs indépendantes pour chaque Yj et tout peut arriver.
 81

Fréquence (IRN)

aft
10,000

5,000

∆
-150 -100 -50 0 50 100 150

Fréquence (CRN)

2,000

1,000

0
Dr
0 5 10 15 20
∆
 82

Un modèle d’inventaire très simple

Xj = niveau d’inventaire au matin du jour j;
Dj = demande durant le jour j, de loi uniforme sur {0, 1, . . . , L};

aft
min(Dj , Xj ) ventes durant le jour j;
Yj = max(0, Xj − Dj ) niveau d’inventaire à la fin du jour j;

Politique (s, S) de gestion d’inventaire: Si Yj < s, on commande S − Yj items.

Chaque commande arrive le lendemain matin avec probabilité p.
Si Yj < s, alors Xj+1 = S avec probabilité p, sinon Xj+1 = Yj .
On pose aussi X0 = S. {Xj , j ≥ 0} est une chaı̂ne de Markov.

Dr
 82

Un modèle d’inventaire très simple

Xj = niveau d’inventaire au matin du jour j;
Dj = demande durant le jour j, de loi uniforme sur {0, 1, . . . , L};

aft
min(Dj , Xj ) ventes durant le jour j;
Yj = max(0, Xj − Dj ) niveau d’inventaire à la fin du jour j;

Politique (s, S) de gestion d’inventaire: Si Yj < s, on commande S − Yj items.

Chaque commande arrive le lendemain matin avec probabilité p.
Si Yj < s, alors Xj+1 = S avec probabilité p, sinon Xj+1 = Yj .
On pose aussi X0 = S. {Xj , j ≥ 0} est une chaı̂ne de Markov.

Dr
[Revenu pour le jour j] = [ventes − coûts d’inventaire − coûts de commande]
= c · min(Dj , Xj ) − h · Yj − (K + k · (S − Yj )) · I[la commande arrive].
On cherche à optimiser les valeurs (s, S). On va les comparer avec les mêmes uniformes.
Deux séquences (“streams”) de nombres aléatoires, une sous-séquence pour chacune des n
répétitions. Mêmes séquences et sous-séquences pour toutes les politiques (s, S).
 83

Modèle d’inventaire: code Java avec la librairie SSJ pour simuler m jours
public double simulateOneRun (int m, int s, int S,

aft
RandomStream streamDemand, RandomStream streamOrder) {
// Simulates inventory model for m days, with the (s,S) policy.
int Xj = S, Yj; // Stock in morning and in evening.
double profit = 0.0; // Cumulated profit.
for (int j = 0; j < m; j++) {
Yj = Xj - [Link] (0, L); // Subtract day demand.
if (Yj < 0) Yj = 0; // Lost demand.
profit += c * (Xj - Yj) - h * Yj;
if ((Yj < s) && ([Link]() < p)) {
// We have a successful order.

}
}
Xj = S;
} else
Xj = Yj; Dr
profit -= K + k * (S - Yj);

return profit / m; // Average profit per day.

 84

Comparaison de p politiques avec CRNs

// Simulate n runs with CRNs for p policies (s[k], S[k]), k=0,...,p-1.

aft
RandomStream streamDemand = new MRG32k3a();
RandomStream streamOrder = new MRG32k3a();
for (int k = 0; k < p; k++) { // Perform n independent runs for policy k.
for (int i = 0; i < n; i++) {
stat_profit[k, i] = simulateOneRun (m, s[k], S[k], streamDemand, streamOrder);
// Advance both streams to next substream, after each run.
[Link]();
[Link]();
}
[Link]();

Dr
[Link]();

// Print and plot results ...

 84

Comparaison de p politiques avec CRNs

// Simulate n runs with CRNs for p policies (s[k], S[k]), k=0,...,p-1.

aft
RandomStream streamDemand = new MRG32k3a();
RandomStream streamOrder = new MRG32k3a();
for (int k = 0; k < p; k++) { // Perform n independent runs for policy k.
for (int i = 0; i < n; i++) {
stat_profit[k, i] = simulateOneRun (m, s[k], S[k], streamDemand, streamOrder);
// Advance both streams to next substream, after each run.
[Link]();
[Link]();
}
[Link]();

Dr
[Link]();

// Print and plot results ...

MRG32k3a avec streams et substreams en Python: [Link]
En MATLAB: [Link]
On pourrait faire exécuter ces pn simulations sur des processeurs parallèles et obtenir exactement les mêmes
résultats, si on utilise les mêmes séquences et sous-séquences.
 85
Comparison with independent random numbers
[Nabil Kemerchou made the plots.]

aft
156 157 158 159 160 161 162 163 164 165 166 167
50 37.94537 37.94888 37.94736 37.95314 37.95718 37.97194 37.95955 37.95281 37.96711 37.95221 37.95325 37.92063
51 37.9574 37.9665 37.95732 37.97337 37.98137 37.94273 37.96965 37.97573 37.95425 37.96074 37.94185 37.93139
52 37.96725 37.96166 37.97192 37.99236 37.98856 37.98708 37.98266 37.94671 37.95961 37.97238 37.95982 37.94465
53 37.97356 37.96999 37.97977 37.97611 37.98929 37.99089 38.00219 37.97693 37.98191 37.97217 37.95713 37.95575
54 37.97593 37.9852 37.99233 38.00043 37.99056 37.9744 37.98008 37.98817 37.98168 37.97703 37.97145 37.96138
55 37.97865 37.9946 37.97297 37.98383 37.99527 38.00068 38.00826 37.99519 37.96897 37.96675 37.9577 37.95672
56 37.97871 37.9867 37.97672 37.9744 37.9955 37.9712 37.96967 37.99717 37.97736 37.97275 37.97968 37.96523
57 37.97414 37.97797 37.98816 37.99192 37.9678 37.98415 37.97774 37.97844 37.99203 37.96531 37.97226 37.93934
58 37.96869 37.97435 37.9625 37.96581 37.97331 37.95655 37.98382 37.97144 37.97409 37.96631 37.96764 37.94759
59 37.95772 37.94725 37.9711 37.97905 37.97504 37.96237 37.98182 37.97656 37.97212 37.96762 37.96429 37.93976
60 37.94434 37.95081 37.94275 37.95515 37.98134 37.95863 37.96581 37.95548 37.96573 37.93949 37.93839 37.9203
61 37.922 37.93006 37.92656 37.93281 37.94999 37.95799 37.96368 37.94849 37.954 37.92439 37.90535 37.93375

IRN

Dr
37.84-37.86
38.02
38
37.98
37.96
37.94
37.92
37.9
37.88
37.86
37.84
156

37.86-37.88
157 158

37.88-37.9
159 160

37.9-37.92
161 162

37.92-37.94
163 164

37.94-37.96
165

37.96-37.98
166
167

37.98-38
50
52
56
54

38-38.02
60
58

38.02-38.02
 86
Comparison with CRNs

aft
156 157 158 159 160 161 162 163 164 165 166 167
50 37.94537 37.94888 37.95166 37.95319 37.95274 37.95318 37.94887 37.94584 37.94361 37.94074 37.93335 37.92832
51 37.9574 37.96169 37.96379 37.96524 37.96546 37.96379 37.96293 37.95726 37.95295 37.94944 37.94536 37.93685
52 37.96725 37.97117 37.97402 37.97476 37.97492 37.97387 37.971 37.96879 37.96184 37.95627 37.95154 37.94626
53 37.97356 37.97852 37.98098 37.98243 37.98187 37.98079 37.97848 37.97436 37.97088 37.96268 37.95589 37.94995
54 37.97593 37.98241 37.98589 37.98692 37.98703 37.98522 37.9829 37.97931 37.97397 37.96925 37.95986 37.95186
55 37.97865 37.98235 37.9874 37.9894 37.98909 37.9879 37.98483 37.98125 37.97641 37.96992 37.96401 37.95343
56 37.97871 37.98269 37.98494 37.98857 37.98917 37.98757 37.98507 37.98073 37.97594 37.96989 37.96227 37.95519
57 37.97414 37.98035 37.98293 37.98377 37.98603 37.98528 37.98239 37.97858 37.97299 37.96703 37.95981 37.95107
58 37.96869 37.97207 37.97825 37.97944 37.97895 37.97987 37.97776 37.97358 37.96848 37.9617 37.95461 37.94622
59 37.95772 37.96302 37.9663 37.97245 37.97234 37.97055 37.9701 37.96664 37.96122 37.95487 37.94695 37.93871
60 37.94434 37.94861 37.95371 37.95691 37.96309 37.96167 37.9586 37.95678 37.95202 37.9454 37.93785 37.92875
61 37.922 37.93169 37.93591 37.94085 37.94401 37.95021 37.94751 37.94312 37.94 37.93398 37.92621 37.91742

CRN

Dr 38
37.98
37.96
37.94
37.92
37.9
37.88
156 157

37.88-37.9
158 159 160

37.9-37.92
161 162

37.92-37.94
163 164

37.94-37.96
165 166

37.96-37.98
167

37.98-38
50
54

52
56
60
58
 87

Estimation de dérivées (sensibilité)

Supposons que la quantité µ que l’on veut estimer dépend d’un paramètre θ dans le modèle,

aft
µ = µ(θ) = Eθ [X (θ)], et on veut estimer
∂µ(θ) ∂Eθ [X (θ)] Eθ1 +δ [X (θ1 + δ)] − Eθ1 [X (θ1 )]
µ′ (θ1 ) = = = lim .
∂θ θ=θ1 ∂θ θ=θ1 δ→0 δ

Ce θ peut apparaitre directement dans le coût X ou encore dans les lois de probabilité.

Dr
 87

Estimation de dérivées (sensibilité)

Supposons que la quantité µ que l’on veut estimer dépend d’un paramètre θ dans le modèle,

aft
µ = µ(θ) = Eθ [X (θ)], et on veut estimer
∂µ(θ) ∂Eθ [X (θ)] Eθ1 +δ [X (θ1 + δ)] − Eθ1 [X (θ1 )]
µ′ (θ1 ) = = = lim .
∂θ θ=θ1 ∂θ θ=θ1 δ→0 δ

Ce θ peut apparaitre directement dans le coût X ou encore dans les lois de probabilité.
Pourquoi estimer les dérivées?
(a) Évaluer l’importance des différents paramètres d’un modèle

Dr
(par ex. pour construire un méta-modèle).
(b) Calculer un Intervalle de confiance qui tient compte de l’erreur d’estimation
des paramètres du modèle.
(c) Évaluer l’effet du changement d’un paramètre de décision.
(d) Un estimateur du gradient est souvent requis dans les algorithmes d’optimisation.
 88

On veut estimer
∂Eθ [X (θ)] Eθ1 +δ [X (θ1 + δ)] − Eθ1 [X (θ1 )]

aft
µ′ (θ1 ) = = lim .
∂θ θ=θ1 δ→0 δ

Différences finies: On choisit δ > 0 très petit.

On simule à θ = θ1 pour obtenir un estimateur X1 = X1 (θ1 ) de µ(θ1 ),
puis on simule à θ = θ2 = θ1 + δ pour obtenir un estimateur X2 de µ(θ2 ),
et on estime la dérivée µ′ (θ1 ) par ∆/δ = (X2 − X1 )/δ.

Dr
 88

On veut estimer
∂Eθ [X (θ)] Eθ1 +δ [X (θ1 + δ)] − Eθ1 [X (θ1 )]

aft
µ′ (θ1 ) = = lim .
∂θ θ=θ1 δ→0 δ

Différences finies: On choisit δ > 0 très petit.

On simule à θ = θ1 pour obtenir un estimateur X1 = X1 (θ1 ) de µ(θ1 ),
puis on simule à θ = θ2 = θ1 + δ pour obtenir un estimateur X2 de µ(θ2 ),
et on estime la dérivée µ′ (θ1 ) par ∆/δ = (X2 − X1 )/δ.
Cet estimateur est biaisé, mais le biais β → 0 quand δ → 0. On a aussi

Var[∆/δ] =

Dr Var(X2 − X1 )
δ
Si X1 et X2 sont indépendants, alors
2
=
Var[X1 ] + Var[X2 ] − 2Cov[X1 , X2 ]

Var[∆/δ] ≈ 2Var[X1 ]/δ 2 = Θ(1/δ 2 ) → ∞ quand δ → 0.

δ2
.
 88

On veut estimer
∂Eθ [X (θ)] Eθ1 +δ [X (θ1 + δ)] − Eθ1 [X (θ1 )]

aft
µ′ (θ1 ) = = lim .
∂θ θ=θ1 δ→0 δ

Différences finies: On choisit δ > 0 très petit.

On simule à θ = θ1 pour obtenir un estimateur X1 = X1 (θ1 ) de µ(θ1 ),
puis on simule à θ = θ2 = θ1 + δ pour obtenir un estimateur X2 de µ(θ2 ),
et on estime la dérivée µ′ (θ1 ) par ∆/δ = (X2 − X1 )/δ.
Cet estimateur est biaisé, mais le biais β → 0 quand δ → 0. On a aussi

Var[∆/δ] =

Dr Var(X2 − X1 )
δ
Si X1 et X2 sont indépendants, alors
2
=
Var[X1 ] + Var[X2 ] − 2Cov[X1 , X2 ]

Var[∆/δ] ≈ 2Var[X1 ]/δ 2 = Θ(1/δ 2 ) → ∞ quand δ → 0.

δ2

Idée: utiliser des valeurs aléatoires communes pour avoir Cov[X1 , X2 ] > 0.
.
 89
Parfois, on peut prendre limδ→0 (X2 − X1 )/δ comme estimateur.

aft
Dr
 89
Parfois, on peut prendre limδ→0 (X2 − X1 )/δ comme estimateur.

Détails: Supposons que X (θ) = f (θ, U) où U ∼ U(0, 1)s et que

aft
∂f (θ, U) f (θ + δ, U) − f (θ, U)
X ′ (θ) = f ′ (θ, U) = = lim
∂θ δ→0 δ
existe avec prob.1 à θ1 . Cette dérivée stochastique f ′ (θ, U) est un estimateur sans biais de
µ′ (θ) ssi  
′ def ∂f (θ, U) ? ∂E[f (θ, U)] def ′
E[f (θ, U)] = E = = µ (θ). (1)
∂θ ∂θ

Dr
 89
Parfois, on peut prendre limδ→0 (X2 − X1 )/δ comme estimateur.

Détails: Supposons que X (θ) = f (θ, U) où U ∼ U(0, 1)s et que

aft
∂f (θ, U) f (θ + δ, U) − f (θ, U)
X ′ (θ) = f ′ (θ, U) = = lim
∂θ δ→0 δ
existe avec prob.1 à θ1 . Cette dérivée stochastique f ′ (θ, U) est un estimateur sans biais de
µ′ (θ) ssi  
′ def ∂f (θ, U) ? ∂E[f (θ, U)] def ′
E[f (θ, U)] = E = = µ (θ). (1)
∂θ ∂θ
Condition suffisante: théorème de convergence dominée (TCD) de Lebesgue.

Dr
S’il existe δ1 > 0 et une v.a. Y tels que

sup
δ∈(0,δ1 ]

et E[Y ] < ∞, alors (1) est valide.

|f (θ + δ, U) − f (θ, U)|
δ
≤Y

Si plusieurs paramètres, le vecteur des dérivées est le gradient stochastique.

(2)
 90

Exemple: réseau d’activités stochastique, E[T ].

5 8 puit

aft
Y5 Y10
Y9 Y12

2 4 7
Y4 Y8
Y1
Y2 Y7 Y11

source 0 1 3 6
Y0 Y3 Y6

−1
Yj = Fj,θ j

Dr
(Uj ). Voir diapo 5 pour les Fj,θj .
Certains Yj suivent une loi exponentielle de moyenne θj :
Yj = Yj (θj ) = −θj ln(1 − Uj ).
Certains Yj suivent une loi normale de moyenne θj et écart-type θj /4:
Yj = Yj (θj ) = θj + (θj /4)Zj = θj + (θj /4)Φ−1 (Uj ).
 91

aft
On veut estimer la dérivée de E[T ] p.r. à chaque θj = µj = E[Yj ].
On considère un θj à la fois. On écrit T = fj (θj , U) où U = (U1 , . . . , U13 ).

Dr
 91

aft
On veut estimer la dérivée de E[T ] p.r. à chaque θj = µj = E[Yj ].
On considère un θj à la fois. On écrit T = fj (θj , U) où U = (U1 , . . . , U13 ).

Différences finies avec CRN: On simule à θj et soit X1 la valeur de T obtenue, puis on simule
à θj + δ et soit X2 la valeur de T obtenue. L’estimateur de dérivée est

X2 − X1 fj (θj + δ, U) − fj (θj , U)
= .
δ δ

Dr
Si on veut estimer tout le gradient (la dérivée pour chaque j), on peut garder le même
X1 = fj (θj , U) pour tous les j, mais on doit simuler un nouveau X2 pour chaque j. Donc
(d + 1)n = 14n simulations au total si on répète n fois.
 92
Dérivée stochastique: L’estimateur de dérivée est

X2 − X1 fj (θj + δ, U) − fj (θj , U)
lim = lim = fj′ (θj , U).

aft
δ→0 δ δ→0 δ
Quand δ est suffisamment petit, le plus long chemin demeure le même. Dans ce cas, X2 − X1
est égal au changement sur Yj si l’arc j est sur le plus long chemin, et X2 − X1 = 0 sinon. Par
conséquent, fj′ (θj , U) = Yj′ (θj ) si l’arc j est sur le plus long chemin, et fj′ (θj , U) = 0 sinon.

Dr
 92
Dérivée stochastique: L’estimateur de dérivée est

X2 − X1 fj (θj + δ, U) − fj (θj , U)
lim = lim = fj′ (θj , U).

aft
δ→0 δ δ→0 δ
Quand δ est suffisamment petit, le plus long chemin demeure le même. Dans ce cas, X2 − X1
est égal au changement sur Yj si l’arc j est sur le plus long chemin, et X2 − X1 = 0 sinon. Par
conséquent, fj′ (θj , U) = Yj′ (θj ) si l’arc j est sur le plus long chemin, et fj′ (θj , U) = 0 sinon.

Si Yj suit la loi exponentielle, alors Yj = Yj (θj ) = −θj ln(1 − Uj ), Yj′ (θj ) = − ln(1 − Uj ) et

fj (θj + δ, U) − fj (θj , U) −δ ln(1 − Uj )

0≤

Dr δ
≤
δ
= − ln(1 − Uj ) = Ej ,

où Ej ∼ Expon(1). Le théorème de convergence dominée s’applique car (2) est vérifiée pour
Y = Ej . Donc E[fj′ (θj , U)] = ∂E[T ]/∂θj (sans biais).
 92
Dérivée stochastique: L’estimateur de dérivée est

X2 − X1 fj (θj + δ, U) − fj (θj , U)
lim = lim = fj′ (θj , U).

aft
δ→0 δ δ→0 δ
Quand δ est suffisamment petit, le plus long chemin demeure le même. Dans ce cas, X2 − X1
est égal au changement sur Yj si l’arc j est sur le plus long chemin, et X2 − X1 = 0 sinon. Par
conséquent, fj′ (θj , U) = Yj′ (θj ) si l’arc j est sur le plus long chemin, et fj′ (θj , U) = 0 sinon.

Si Yj suit la loi exponentielle, alors Yj = Yj (θj ) = −θj ln(1 − Uj ), Yj′ (θj ) = − ln(1 − Uj ) et

fj (θj + δ, U) − fj (θj , U) −δ ln(1 − Uj )

0≤

Dr δ
≤

Y = Ej . Donc E[fj′ (θj , U)] = ∂E[T ]/∂θj (sans biais).

δ
= − ln(1 − Uj ) = Ej ,

où Ej ∼ Expon(1). Le théorème de convergence dominée s’applique car (2) est vérifiée pour

Si Yj suit la loi normale, alors Yj = Yj (θj ) = θj + (θj /4)Φ−1 (Uj ) et Yj′ (θj ) = 1 + Φ−1 (Uj )/4.
Sans biais aussi (facile à vérifier).
 93

Exemple: réseau d’activités stochastique, P[T > x].

On veut maintenant estimer la dérivée de P[T > x] p.r. à θj .
L’estimateur de P[T > x] est fj (θj , U) = I[T > x]. Ne peut prendre que les valeurs 0 et 1.

aft
Dr
 93

Exemple: réseau d’activités stochastique, P[T > x].

On veut maintenant estimer la dérivée de P[T > x] p.r. à θj .
L’estimateur de P[T > x] est fj (θj , U) = I[T > x]. Ne peut prendre que les valeurs 0 et 1.

aft
La dérivée fj′ (θj , U) est toujours soit 0, soit pas définie (survient avec prob. 0).
On a donc P[fj′ (θj , U) = 0] = 1. Estimateur biaisé de µ′ (θj ) = ∂P[T > x]/∂θj .
Ici le théorème de convergence dominée ne s’applique pas car
supδ>0 [fj (θj + δ, U) − fj (θj , U)]/δ n’est pas intégrable.
Le problème vient du fait que fj est discontinue en θj .

Dr
 93

Exemple: réseau d’activités stochastique, P[T > x].

On veut maintenant estimer la dérivée de P[T > x] p.r. à θj .
L’estimateur de P[T > x] est fj (θj , U) = I[T > x]. Ne peut prendre que les valeurs 0 et 1.

aft
La dérivée fj′ (θj , U) est toujours soit 0, soit pas définie (survient avec prob. 0).
On a donc P[fj′ (θj , U) = 0] = 1. Estimateur biaisé de µ′ (θj ) = ∂P[T > x]/∂θj .
Ici le théorème de convergence dominée ne s’applique pas car
supδ>0 [fj (θj + δ, U) − fj (θj , U)]/δ n’est pas intégrable.
Le problème vient du fait que fj est discontinue en θj .

Dr
Solution: On peut régler ce problème en utilisant CMC:
On remplaçe l’indicateur I[T > x] par la probabilité conditionnelle

P[T > x | {Yj , j ̸∈ L}] = P[T > x | Y0 , Y1 , Y2 , Y3 , Y7 , Y10 , Y11 , Y12 ],

qui est continue en θj . Voir diapos 76–77.

La dérivée stochastique de cet estimateur donne un estimateur sans biais de la dérivée.