Statistiques Bayésiennes

Hicham Janati
hjanati@[Link]
1. Pourquoi Monte-Carlo ? (Exemple de modèle hiérarchique)
2. Introduction à la méthode Monte-Carlo (historique, PRNG)

3. Algorithmes de simulation i.i.d (PRNG, transformation, rejet)

4. Méthodes MCMC (Gibbs, Metropolis)

5. Diagonstics de convergence MCMC

6. Méthodes MCMC avancées (Langevin, HMC, NUTS)

2
 Pourquoi les
II. Méthodes de Monte-Carlo
?

3
Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique

On souhaite modéliser la fréquence des sinistres d’un ensemble de conducteurs:

TD1

Modèle (vraisemblance) a priori

i ⇠ Gamma(↵, )
Nombre de sinistres par années N i | i ⇠ P( i )
issu de données historiques
sur tous les clients

Or les conducteurs peuvent être bons ou mauvais:

Inconvénients: tous les conducteurs ont la même loi a priori

Sur-estimer le risque des bons conducteurs / sous-estimer le risque des mauvais conducteurs
On note z i = 1 (bon) et z i = 0 (mauvais)

4
Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique

On souhaite modéliser la fréquence des sinistres d’un ensemble de conducteurs:

TD1

Modèle (vraisemblance) a priori

i ⇠ Gamma(↵, )
Nombre de sinistres par années N i | i ⇠ P( i )
issu de données historiques
sur tous les clients

Or les conducteurs peuvent être bons ou mauvais:

Inconvénients: tous les conducteurs ont la même loi a priori

Sur-estimer le risque des bons conducteurs / sous-estimer le risque des mauvais conducteurs
On note z i = 1 (bon) et z i = 0 (mauvais)
Mais les z i ne sont pas observées...
Comment peut-on adapter le modèle ?
5
 Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique
N i| i ⇠ P( i )
On note z i = 1 (bon) et z i = 0 (mauvais)
Comment peut-on modéliser la variable z i ?
On modélise la variable z i |p ⇠ Bernoulli(p)
Lois a priori des bons et des mauvais di érentes:
groupe a une loi a priori: i |z i = 0 ⇠ Gamma(↵0 , 0) et i |z i = 1 ⇠ Gamma(↵1 , 1)

Comment choisir les paramètres ?

↵ 0 , ↵1 , 0, 1 ⇠ Gamma(⌧, )
L’hyperprior Gamma “lie”
Avec des hyperpriors: les deux groupes
p ⇠ (a, b)

Comment choisir les hyper-paramètres ?

Ceci est un exemple de Modèle bayésien hiérarchique

6
ff
 Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique

Modèle bayésien hiérarchique comme un PGM (Probabilistic graphical model)

On représente les dépendances entre variables aléatoires avec une èche:

N i| i ⇠Données
P( i ) de l’individu
N i| i ⇠ P( i )
Paramètres individuels
N i | i ⇠ P( i )
Et chaque groupe a une loi a priori: i |z i = 0 ⇠ Gamma(↵0 , 0 ) et i |z i = 1 ⇠ Gamma(↵1 , 1 )
priori: i |z i = 0 ⇠ Gamma(↵0 , 0 ) et i |z i = 1 ⇠ Gamma(↵1 , 1 )
Et chaque groupe a une loi a priori: i |z i = 0 ⇠ Gamma(↵0 , 0 ) et i |z i = 1 ⇠ Gamma(↵1 , 1 )
pe
upe Variables
aEtune
aEt
chaque
une
chaque
loiloi latentes
a groupe
priori:
a groupea(non-observées)
priori: |z
iune |z=iloi
a iiune = ⇠
0loi a⇠
a0 priori:
Gamma(↵
priori:
Gamma(↵ i |zi0i|z
,= 0) 0⇠
0i ,0= 0) ⇠
et |zii|z=i =
Gamma(↵
etiGamma(↵ 1⇠01, ⇠
Gamma(↵ ) eti |zi1i|z
0 ,0 Gamma(↵
) 0et ,= 1) 1⇠
1i ,1= 1) ⇠Paramètres
Gamma(↵
Gamma(↵ 1 ,collectifs
1 ,1 ) 1 )
élise la variable z i |p ⇠ Bernoulli(p)
n modélise la variable z i |p ⇠ Bernoulli(p) ↵0 , ↵1 , 0 , 1 ⇠ Gamma(⌧, )
↵ , ↵ , ↵, , ↵⇠, Gamma(⌧,
0 1 00 1 1 0 , ⇠ Gamma(⌧,
1 ) )
p ⇠ (a, b) Hyper-paramètres
p ⇠ p(a,
⇠ b)(a, b)

7
fl
 Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique

Modèle bayésien hiérarchique comme un PGM (Probabilistic graphical model)

On représente les dépendances entre variables aléatoires avec une èche:

N i | i ⇠ P( i )
Ne pas oublier que chaque
N i| i ⇠…P( i )N i | i …
⇠ P( N
i ) i | i ⇠ P( i )
N i | i ⇠ P( i ) individu a son paramètre
individuel…
… … )
N i | i ⇠ P( i )N i | i ⇠ P( Et chaque individu peut avoir
i i | i ⇠ P( i )
N
plusieurs observations …

iroupe
Et achaquea une
priori: loii =
groupe
i |z … Gamma(↵
a apriori:
0 une
⇠ loi ia|zpriori:
i =… 00, ⇠0 )Gamma(↵
i |z
eti =i |z 0 i⇠=0Gamma(↵ et i |z
,1 0⇠) Gamma(↵ 0 ,i = , ⇠1 )Gamma(↵
0 )11et i |z i = 1 ⇠1Gamma(↵, 1) 1, 1)
peupe
aEtune
aEt
chaque
une
chaque
loiloi
a groupe
priori:
a groupe
priori: |z
a iune |z=iloi
a iiune = ⇠
0loi a⇠
a0 priori:
Gamma(↵
priori:
Gamma(↵ i |zi0i|z
,= 0) 0⇠
0i ,0= 0) ⇠
et |zii|z=i =
Gamma(↵
etiGamma(↵ 1⇠ 01, ⇠
Gamma(↵ ) eti |zi1i|z
0 ,0 Gamma(↵
) 0et ,= 1) 1⇠
1i ,1= 1) ⇠
Gamma(↵
Gamma(↵
1 , 1 ,1 ) 1 )

n modélise la variable z i |p ⇠ Bernoulli(p)

↵0 , ↵1 , 0↵,0 , 1↵1⇠, Gamma(⌧,
0 , 1 ⇠ Gamma(⌧,
) )

p ⇠ p(a,
⇠ b)(a, b) Aussi appelé un “Réseau bayésien”

8
fl
 Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique

Modèle bayésien hiérarchique comme un PGM (Probabilistic graphical model)

On représente les dépendances entre variables aléatoires avec une èche:

N i| i ⇠Données
P( i ) de l’individu
N i| i ⇠ P( i )
Paramètres individuels
N i | i ⇠ P( i )
Et chaque groupe a une loi a priori: i |z i = 0 ⇠ Gamma(↵0 , 0 ) et i |z i = 1 ⇠ Gamma(↵1 , 1 )
priori: i |z i = 0 ⇠ Gamma(↵0 , 0 ) et i |z i = 1 ⇠ Gamma(↵1 , 1 )
Et chaque groupe a une loi a priori: i |z i = 0 ⇠ Gamma(↵0 , 0 ) et i |z i = 1 ⇠ Gamma(↵1 , 1 )
pe
upe Variables
aEtune
aEt
chaque
une
chaque
loiloi latentes
a groupe
priori:
a groupea(non-observées)
priori: |z
iune |z=iloi
a iiune = ⇠
0loi a⇠
a0 priori:
Gamma(↵
priori:
Gamma(↵ i |zi0i|z
,= 0) 0⇠
0i ,0= 0) ⇠
et |zii|z=i =
Gamma(↵
etiGamma(↵ 1⇠01, ⇠
Gamma(↵ ) eti |zi1i|z
0 ,0 Gamma(↵
) 0et ,= 1) 1⇠
1i ,1= 1) ⇠Paramètres
Gamma(↵
Gamma(↵ 1 ,collectifs
1 ,1 ) 1 )
élise la variable z i |p ⇠ Bernoulli(p)
n modélise la variable z i |p ⇠ Bernoulli(p) ↵0 , ↵1 , 0 , 1 ⇠ Gamma(⌧, )
↵ , ↵ , ↵, , ↵⇠, Gamma(⌧,
0 1 00 1 1 0 , ⇠ Gamma(⌧,
1 ) )
p ⇠ (a, b) Hyper-paramètres (non aléatoires)
p ⇠ p(a,
⇠ b)(a, b)

n n n
Question: trouver la loi a posteriori jointe de tous les paramètres ({ i }i=1 , {z i }i=1 , p, ↵0 , ↵1 , 0, 1 ) |{N i }i=1 .

9
fl
 Pourquoi Monte-Carlo ? Modèle bayésien hiérarchique

La loi jointe ({ i }ni=1 , {z i }ni=1 , p, ↵0 , ↵1 , 0, 1 ) |{N }n

i i=1 est intractable !

( )
⌘ Yn ⇣ ⌘ ⇣ ⌘ ⇣ ⌘
, 1 , {zi }, { i } | {Ni } / P Ni | i P i | zi , ↵ 0 , 0 , ↵ 1 , 1 P z i | p P p P ↵ 0 P 0 P ↵1 P 1 .
i=1

( ✓ ◆1 ) "" #
n
Y Ni ↵0 zi ✓ ↵1 ◆z i
i e i
⌧ ⌧ 11 e ↵00 ⌧ 11 0 0 0
↵10 ↵1⌧ 0 ⌧10 1 0 0
/ 0
i
↵0 1
e 0 i 1
i
↵1 1
e 1 i zi
p (1 p) 1 zi
p a 1
(1 p) b 1 ↵
↵00 e ↵
00 ee ↵1 1 1 e1
00⌧⌧
↵ e 1 1 e e 1
.1
i=1
N i! (↵0 ) (↵1 )

n n
Comment obtenir les lois a posteriori d’intérêt: j |{N i }i=1 ? z j |{N i }i=1 ?

Il faut marginaliser: intégrer la loi jointe par rapport à

toutes les autres variables: impossible analytiquement !

Solution: Obtenir des échantillons “simulés” issus de la loi a posteriori et construire des estimateurs empiriques

Et pour cela on utilise les méthodes de simulation “Monte-Carlo”

10
 Pourquoi les
2. Méthodes de Monte-Carlo
?

1. Modèles Bayésiens: Simuler à partir d’une loi a posteriori intractable

2. Résoudre des équations di érentielles / intégrales en Physique / Finance / Météo
3. Incontournables en IA / ML (Bayesian Optimization, Generative AI, Reinforcement
learning …)

11
ff
II Méthodes de Monte-Carlo
1. Introduction

12
 Calcul d’intégrale pre-WWII
Nous sommes en 1902
?
Cible X

? ?
Z = XZ =XX X
Vous
On note d = 30 le diamètre de la cib
?
On note X et X les coordonnées du projectile et de la cible.
? 2
On modélise X ⇠ N (X , ) avec = 50.
On note d = 30 le diamètre de la cible.

Pour détruire la cible, il faut en moins 40 impacts directs, combien de projectiles devez-vous lancer ?

def X X ?
On pose Z = ⇠ N (0, 1). P(Impact)
P(Impact)==P(|X
P(|X ??
X | d)==P(|Z|
X |d) dd
P(|Z| ))== (( )) dd
(( dd
))
d 30
\ = =2 ( ) 1=2 ( 50 ) 1 ⇡ 0.45
P(Impact) 40
N )N = \
40
⇡ 40
0.45 ⇡ 89
P(Impact) Comment on calcule 30
( 50 ) en 1902 ?
13
 Calcul d’intégrale pre-WWII
30
Comment on calcule ( 50 ) en 1902 ?
Z x
def 1 x2
(x) = P(Z  x) = p e 2 2 dx
2⇡ 1

Comment calculer ⇡ ? La racine ? l’exponentielle ? l’intégrale ?

x3 x5 x7
1. ⇡ = arctan(1)[Link]⇡ arctan(x)
= arctan(1)= xet arctan(x)
+ = x + ·
x ·
3
· + x 5
x 7
+ · · ·
x33
x55
x77
3 5 7
1. ⇡ = arctan(1) et arctan(x) = x⇣ 3 + ⌘5 + · · ·
ler ⇡ ?[Link] racine ? l’exponentielle ? l’intégrale ?
1 ⇣ 3 A ⌘ 5
7⇣
7 ⌘p p
La=suite de Newton a = a x+ x x
converge
1 + · ·vers A pA.
1. ⇡ arctan(1) et arctan(x) =
2. La suite de Newton
n+1 2
1
x n +
3 aan+1
An 5 = 2 7 an + an
· converge vers A.
2. La suite de Newton an+1 = 2 ⇣an x+3 an ⌘ 5 converge
7 vers A.
1. ⇡ = arctan(1) et arctan(x) =1 x + x
Px
+ · · · p
2.
3. La suite de Newton
Approcher l’intégralea = a
avec l’intégrale
des 3
+ A 5 7
converge
⇡des
10 vers A. xP2
2
2 ⇣rectangles exp( ⇡2x22 )
n+1 n i 10 x
tielle ?3.l’intégrale ? 3. Approcher aavec
n⌘ Pi=1 rectangles
10 p i=1 exp( 2
i
2)
Approcher l’intégrale avec des 1 rectangles A ⇡ i=1 exp( 2
i
2)
2. La suite de Newton an+1 = 2 an +x2an xconverge 3 P vers A. 2
4.
3. Calculer les exp: exp(x)
Approcher l’intégrale
4. Calculer =
avec 1 +
lesdes
exp: x +
rectangles
exp(x)
2!
2
+ = ⇡
3!
3 1++ · x·
10 · + x 2
exp(
+ x 3 x
+ i
2·) · ·
i=1 2
4. Calculer les exp: exp(x) = 1 + x + x2! + x3! + P10· · · 2! 3!
2
x
3. Approcher l’intégrale avec des rectangles x⇡ exp( 2 2 )
2 3 i
x
Calculer les exp: exp(x) = 1 + x + 2! + 3! + · · ·
acine ? 4.l’exponentielle ? l’intégrale ? i=1

x2 x3 x
4. Calculer les exp: exp(x) = 1 + x + 2! + 3! + · · ·
Beaucoup de calculs à la main…
14
Calcul d’intégrale pre-WWII
“Un computer” à l’époque était un travail:

1] produire des livres avec

des tables de calcul
précises cosinus, sinus,
exp, logarithme, racines,
puissances, les constantes,
intégrales…

2] Résoudre manuellement (algorithmes) des équations di érentielles en

mécanique (trajectoires de projectiles) en utilisant des calculatrices
mécaniques de bureau

3] Produire des tables de calculs des trajectoires en fonction de

l’angles de tir, vitesse et direction du vent …

15
ff
ENIAC et méthodes de quadrature

1946: Premier ordinateur électronique ENIAC

(Electronic Numerical Integrator and Computer)

“Human computers” ENIAC

Calcul d’une trajectoire: 40 heures 30 secondes

Et si on utilisait ENIAC pour calculer des intégrales doubles ? triples ?

Applications: Prévision météo, Finance, Physique nucléaire …

La première prévision météo (ENIAC) du lendemain a nécessité 24h de calcul !

Les méthodes de quadrature sont lentes en grande dimension

16
L’idée “Monte-Carlo”

Stan Ulam: Physicien / mathématicien membre du Manhattan project (bombe atomique)

“Au lieu d’étudier tous les cas possibles, je peux jouer

1000 parties et calculer le pourcentage de parties
gagnées…”

“Et si ENIAC peut faire ça à ma place …”

ATTENDS ! Au lieu de calculer toutes les trajectoires de

neutrons possibles avec les quadratures, je pourrais les
(r, ⌦, E) = q(r, ⌦, E) simuler et avoir une approximation de:
À l’hôpital, essaie de calculer la
probabilité de gagner au solitaire en Z 1 Z XL ⇣ XM ⌘
2` + 1
+ A`,m (r) exp B`,m (r) |E 0 E| P` ⌦0 ·⌦ (r, ⌦0 , E 0 ) d⌦0 dE 0 .
analysant toutes les combinaisons de 0 4⇡ 4⇡ m=1
`=0
cartes possibles
L’ordinateur était toujours considéré comme un “calculateur déterministe”. L’utiliser pour
“simuler” des nombres selon des “probabilités” était une idée révolutionnaire
17
 L’idée “Monte-Carlo”
Z 1
Ulam, Metropolis et Von Neumann proposent une méthode g(x)dx = EX⇠U ([0,1]) (g(X
de calcul d’une intégrale quelconque de la forme: 0

Comment peut-on écrire cette intégrale comme une espérance ?

Z 1 X N
1
g(x)dx = EX⇠U ([0,1]) (g(X)) ⇡ Comment g(X i)
approximer cette espérance ?
0 N i=1
N
1 X
x = EX⇠U ([0,1]) (g(X)) ⇡ g(X i ) Avec X 1 , . . . , X N i.i.d ⇠ U([0, 1])
N i=1

Metropolis propose de nommer leur article “The Monte-Carlo Method” inspiré par
l’oncle d’Ulam qui était accro au célèbre casino de Monte-Carlo.

Comment générer X 1 , . . . , X N i.i.d ⇠ U([0, 1]) ?

18
Le premier PRNG

Comment générer x ⇠ U([0, 1]) ?

Première idée de Newman:

1. Choisir un nombre à 10 chi res: 4924748149

2. Calculer son carré: 24253144331078926201

3. Prendre 10 chi res au milieu: 24253144331078926201

4. Diviser par 1010: x = 1443310789 / 10000000000

5. Pour un second chi re, calculer 14433107892 … et ainsi de suite

La suite à 10 chi res est une suite périodique (déterministe) avec un comportement qui “semble” aléatoire

On appelle ce type d’algorithme: PRNG (Pseudo-random numbers generators)

19
ff
ff
ff
ff
Les PRNG aujourd’hui

1. Les
Se PRNG
basentd’aujourd’hui:
sur un mélange F d’opérations (addition, multiplication, reste de division, shift de bits . . . )
1.
1. Se
Se basent
basent sur
sur un
un mélange
mélange F
F d’opérations
d’opérations (addition,
(addition, multiplication,
multiplication, reste
reste de
de division,
division, shift
shift de
de bits
bits .. .. .. ))
1. Se
1. Se basent
basent sur sur unun mélange
mélange F F d’opérations
d’opérations (addition,
(addition, multiplication,
multiplication, reste reste de de division,
division, shift
shift dede bits
bits .. .. .. ))
2. Définissent toujours une suite d’entiers périodique xn+1 = F (xn ) avec xi 2 [0, M ], avec M très grand et une pério
2.
2. Définissent
Définissent toujours
toujours une
une suite
suite d’entiers
d’entiers
12 périodique
périodique x
x n+1
n+1 =
= F
F (x
(x n
n )) 18
avec
avec x
x ii 2
2 [0,
[0, M
M ],
], avec
avec M
M très
très grand
grand et
et une
une période
période
2. très très grande
2. Définissent
Définissent toujours
toujours (même
une
une avec d’entiers
suite
suite 10
12 chi↵res
d’entiers / seconde,
périodique
périodique x
x n+1 il=
=faut
F
F (x
(x 10)) avec
18
n années
avec x
x i 2
2 pour
[0,
[0, M
M voir
],
], avec
avectoute
M
M la séquence).
très
très grand
grand et
et une
une période
période
très
très très
très grande
grande (même
(même avec
avec 10
10 12
12 chi↵res
chi↵res /
/ seconde,
seconde, il
n+1
il faut
faut 10
10 n années
18
18 années ipour
pour voir
voir toute
toute la
la séquence).
séquence).
très très
très très grande
grande (même avec
xi (même avec 10
10 chi↵res
12 chi↵res // seconde,
seconde, il il faut
faut 10
10 années
18 années pour
pour voir
voir toute
toute la
la séquence).
séquence).
3. Définir ui = xxM permet d’obtenir des nombres qui “semblent” suivre
i permet d’obtenir des nombres qui “semblent” suivre une loi uniforme U ([0, 1]).
une loi uniforme U ([0, 1]).
3.
3. Définir
Définir u
u ii =
= x
M
i
i
i permet d’obtenir des nombres qui “semblent” suivre une loi uniforme U ([0, 1]).
3. Définir
3. Définir u uii == Mi permet
M
x permet d’obtenir
d’obtenir desdes nombres
nombres qui qui “semblent”
“semblent” suivre suivre une une loiloi uniforme
uniforme U U ([0,
([0, 1]).
1]).
4.
4. Il
Il existe
existe plusieurs
M
plusieurs testsqui
tests quipermettent
permettent de de mesurer
mesurer la la qualité
qualité statistique
statistique d’un d’un
PRNGPRNG (corrélation,
(corrélation, biais biais
de de sélection.
sélection...)
4. Il
4. Il existe
Il existe plusieurs
existe plusieurs
plusieurs teststests qui
tests qui permettent
qui permettent
permettent de de mesurer
de mesurer
mesurer la la qualité
la qualité statistique
qualité statistique
statistique d’un d’un PRNG
d’un PRNG (corrélation,
PRNG (corrélation, biais
(corrélation, biais
biais de de sélection...)
de sélection...)
sélection...)
4.
5. On
5. On appelle
appellexx0lalagraine graine(seed)
(seed)duduPRNG:
PRNG: si si
on on
la la connait,
connait, on on
peut peut déterminer
déterminer toutetoute
la la suite.
suite.
5. On
5. On appelle
On appelle
appelle x x 0
0
x00 lala graine
la graine (seed)
graine (seed) du
(seed) du PRNG:
du PRNG:
PRNG: si si on
si on la
on la connait,
la connait,
connait, on on peut
on peut déterminer
peut déterminer toute
déterminer toute la
toute la suite.
la suite.
suite.
5.
6.
6.
6. Fixer
Fixer
Fixer la
la
la seed
seed
seed permet
permet
permet de
dede garantir
garantir
garantir la
lala reproducibilité
reproducibilité
reproducibilité d’un
d’und’un programme
programme
programme avec
avec avec
des
des des nombres
nombres
nombres (pseudo)
(pseudo)
(pseudo) aléatoires.
aléatoires.
aléatoires.
6. Fixer
6. Fixer la la seed
seed permet
permet de de garantir
garantir lala reproducibilité
reproducibilité d’un d’un programme
programme avec avec desdes nombres
nombres (pseudo)
(pseudo) aléatoires.
aléatoires.
7.
7.
7. Les
Les
Les PRNG
PRNG
PRNG les
les
les plus
plus
plus connus
connus
connus sont:
sont:
sont:
7. Les
7. Les PRNG
PRNG les les plus
plus connus
connus sont:
sont:
(a) Mersenne
(a)
(a) MersenneTwister
Mersenne Twister(1997):
Twister (1997):utilisé
(1997): utilisépar
utilisé par par défaut
défaut
défaut Numpy
Numpy
Numpy <=
<= <=
1.16.
1.16. 1.16 [Link]
(a) Mersenne
(a) Mersenne Twister Twister (1997):
(1997): utilisé
utilisé par
par défaut
défaut Numpy
Numpy <= <= 1.16.
1.16.
(b) PCG64
(b)
(b) PCG64(2014):
PCG64 (2014):adopté
(2014): adoptépar
adopté parNumpy
par Numpy>=
Numpy >=>=1.171.17
1.17 en
en en 2019. [Link].default_rng
2019.
2019.
(b) PCG64
(b) PCG64 (2014):
(2014): adopté
adopté par
par Numpy
Numpy >= >= 1.17
1.17 enen 2019.
2019.
8. Ne doivent pas être utilisés en cryptographie (générer des clés privées / mots de passe)

20
 Algorithmes de simulation

Les PRNG permettent de générer des échantillons i.i.d X 1 , . . . , X n ⇠ U([0, 1]).

Comment peut-on générer les autres lois ?

U ([a,? b]) ? n)? 2

Qu’en est-il de X 1 , . . . , X n ⇠ B(p) Bin(p, (↵, ) ? P( ) ? N (µ, ) ? ...
2
Qu’en est-il de X 1 , . . . , X n ⇠ B(p) ? Bin(p, n)? ? P( ) ?théorie
(↵, )Résultats N (µ, des) ? ...
Qu’en est-il
Qu’endeest-il . . ,XX1n, .⇠
X 1 , .de . . B(p)
, Xn ⇠ B(p) ? n)?
? Bin(p, Bin(p,
(↵,n)?) ?(↵, probabilités
P( )) ?? P(
N (µ,
)? N2 (au tableau
2 / Python)
) ?(µ,... ) ? ...
Basés sur des Exp(? )Bin(p,
Qu’en est-il de X 1 , . . . , X n ⇠ B(p) ? n)? (↵, ) ? P( ) ? N (µ, 2
) ? ...
échantillons i.i.d
Gamma(k,
Qu’en est-il de X 1 , . . . , X n ⇠ B(p) ? Bin(p, k 2 N(↵,
✓) n)? ? ? ) ? P( ) ? N (µ,
2
) ? ...
2 2
⇠ B(p) ?Qu’en
Bin(p,
est-il X 1 ,). .?. ,P(
n)? de(↵, X n )⇠? B(p)
N (µ,? Bin(p,Box-Muller,
) ? ... ) ? P( ) Ziggurat
n)? (↵, Marsaglia, ? N (µ, ) ? ...
f 2
loi Qu’en est-il
quelconque à densité 1 , . . . , X n ⇠=B(p)
de Xnormalisée R
f
Rejection
? Bin(p,
non n)?sampling
(↵, /) g? ?P( ) ? N (µ,
normalisée ) ? ...
Basés sur des f
tééchantillons
normalisée =
Qu’en non R ? non normalisée / g ? Algorithmes MCMC (Markov-Chain
f de X 1 , . . . , X n ⇠ B(p) ? Bin(p, n)? (↵, ) ? P( ) ? N (µ,
est-il
i.i.d Monte-Carlo)
2
) ? ...

21
Algorithmes de simulation Rejection sampling

Les PRNG permettent de générer des échantillons i.i.d X 1 , . . . , X n ⇠ U([0, 1]).

Comment peut-on générer les autres lois ?

Comment simuler des points à

l’intérieur de la surface bleue ?

22
Algorithmes de simulation Rejection sampling

Les PRNG permettent de générer des échantillons i.i.d X 1 , . . . , X n ⇠ U([0, 1]).

Comment peut-on générer les autres lois ?

On simule des échantillons

uniformes dans le carré [0, 1]2

23
Algorithmes de simulation Rejection sampling

Les PRNG permettent de générer des échantillons i.i.d X 1 , . . . , X n ⇠ U([0, 1]).

Comment peut-on générer les autres lois ?

On “rejette” les échantillons à

l’extérieur de la surface.

24
Algorithmes de simulation Rejection sampling

On souhaite adapter cette méthode pour avoir des échantillons qui suivent une densité f: dé nie sur [-3, 3]:

(x, y)

(x, y)

x à garder x à rejeter

Comment peut-on simuler des échantillons à partir de sa distribution avec des échantillons uniformes ?

Simuler uniformément sur tout le rectangle un point (x, y)

Garder x si le point est sous la densité: si y < f(x)

25

fi
Algorithmes de simulation Rejection sampling

(animated)

26
Algorithmes de simulation Rejection sampling

Et si on avait cette densité comme target ?

Ou celle-ci ?

27
Algorithmes de simulation Rejection sampling
On rejette beaucoup trop d’échantillons.
Comment y remédier ?

(animated)

Idée: prendre une courbe bleue plus proche de la densité cible f

28
 Algorithmes de simulation Rejection sampling

Supposons qu’on sait générer à partir d’une loi à densité g:

Rejection sampling
1.
1. Sample x⇠
Sample x ⇠ gg
2. Sample yy ⇠
2. Sample ⇠ U([0,
U([0, M
M g(x)])
g(x)])
On cherche alors une constante M telle que f < Mg: 3. If yy 
3. If  ff (x):
(x): Accept
Accept x
x else
else Reject
Reject x.
x.

0, M g(x)])
ple x ⇠ g
cept x else Reject x.
ple y ⇠ U([0, M g(x)])
 f (x): Accept x else Reject x.
1. Sample x ⇠ g
f: Target distribution 29
g: Proposal distribution
2. Sample y ⇠ U([0, M g(x)])
Algorithmes de simulation Rejection sampling

30
Algorithmes de simulation Rejection sampling

Algorithme de rejet:
1.
1. On
On souhaite
souhaite simuler
simuler à
à partir
partir d’une
d’une densité
densité ff ..
1. On souhaite simuler à partir d’une densité f .
2.
2. On
On cherche
cherche une
une densité
densité gg selon
selon à
à laquelle
laquelle on
on sait
sait générer
générer des
des échantillons.
échantillons.
2. On cherche une densité g selon à laquelle on sait générer des échantillons.
3.
3. Il
Il faut
faut dominer
dominer ff par
par gg en
en trouvant
trouvant M
M >
> 0
0 tel
tel que:
que: ff 
 M
M gg partout.
partout.
3. Il faut dominer f par g en trouvant M > 0 tel que: f  M g partout.
1
4.
4. La
La probabilité
probabilité d’acceptation
d’acceptation est
est donnée
donnée par
par M
1 : on veut le plus petit M possible
1 : on veut le plus petit M possible (borne
(borne serrée)
serrée)
4. La probabilité d’acceptation est donnée par M M : on veut le plus petit M possible (borne serrée)
5.
5. On
On itère
itère par
par la
la suite
suite l’algorithme:
l’algorithme:
5. On itère par la suite l’algorithme:
1. Sample x ⇠ g Quelles sont les limites de cette méthode ?
2. Sample y ⇠ U([0, M g(x)])
1. Il faut connaître la densité f normalisée !
3. If y  f (x): Accept x else Reject x.
2. En grande dimension, la zone de rejet explose

Solution: Les méthodes MCMC

31
II Méthodes de Monte-Carlo
1. Introduction
2. Markov Chain Monte-Carlo (MCMC)
3. Algorithmes MCMC avancés

32
Markov Chains Dé nition
Dé nition: Chaîne de Markov

Soit (X n )n 0 une suite de variables aléatoires à valeurs dans S.

On dit que (X n )n 0 vérifie la propriété de Markov si son état futur n + 1
ne dépend que du présent n et non pas du passé (n 1 .. 0):
P(X
P(Xn+1 = xx || X
n+1 = Xnn =
= xxnn,,......,,X
X00 = = P(X
= xx00)) = P(Xn+1 = xx || X
n+1 = Xnn =
= xxnn))

Quel graphe probabiliste permet de visualiser cette suite ?

P(X n+1==xx| P(X

P(Xn+1 |XXnn=
P(X =x
n+1x
n+1n=,
n=.
, .
x .
X
x |,
.|X
,
n X
X0
n
1n0==xx0
n )
,
0
n =
.
)
, .
= .P(X
,
. P(X
X
, X0 =
n+1
0 =x
n+1=
x0=)
0 x
)=x
=| P(X
X
| P(X
Xn n==x
n+1x
n+1n=)
=
n )xx| X
| Xnn

Dé nition: Chaîne de Markov homogène

On dit que (X n )n est homogène si les variables X n+1 | X n ont la même loi 8n.

Dans toute la suite, les chaînes de Markov sont supposées homogènes

33
fi
fi
fi
Markov Chains Exemple
Soit (X n )n 0 une chaı̂ne de Markov à valeurs dans S = {A, B, C}.

X
X nn peut
peut valoir
valoir A,
A, B
B ou
ou CC et
et X
X n+1
n+1
ne
ne dépend
dépend que
que de
de X
Xnn
avec
avec les
les probabilités
probabilités conditionelles
conditionelles de
de transition:
transition:

To A To B To C
From A 0.5 0.3 0.2
From B 0.3 0.4 0.3
From C 0.2 0.3 0.5

La matrice de transition P

⇥ ⇤>
On note la distribution de X n par le vecteur de probabilités ⇡n = P(X n = A) P(X n = B) P(X n = C) 2 R3 .
1. Déterminez ⇡n+1 en fonction de P et de ⇡n .

2. On suppose que X n a une distribution asymptotique ⇡, comment peut-on l’obtenir ?

34
Markov Chains Stationnarité: cas discret

Soit (X n )n 0 une chaı̂ne de Markov à valeurs dans un ensemble fini S.

Alors:

1. (X nn )nn 0
0 est dite irréductible si 8i, j 2 S, 9n 2 N tel que on peut atteindre j à partir de i après n étapes.
2. (X nn )nn 0
0 est dite apériodique si la chaı̂ne ne présente pas de comportement cyclique.
3. (X nn )nn 0
0 est dite récurrente
ergodique sipositive
elle est àsilalefois
temps de retour àetn’importe
irréductible quel i 2 S est fini.
apériodique.
4. (X
Unenn )chaı̂ne
n
n 00 est(X
dite
n )nergodique si elle
0 ergodique est àtoujours
admet la fois irréductible, apériodique
une distribution stationnaireet⇡récurrente
= limn!+1positive.
⇡n .
5. Une chaı̂ne (X nn )nn 0
0 ergodique admet toujours une distribution stationnaire unique ⇡ = limn!+1
n!+1 ⇡n
n.

Histogramme

des états A, B, C
35 (animated)
Markov Chains Stationnarité: cas continu
Soit (X n )n 0 une chaı̂ne de Markov à valeurs dans Rd .
La matrice de transition P devient un “noyau”: (x, y) 7! K(x, y) qui correspond à la densité de X n+1 = y|X n = x
R R
n ⇠
Si XSi X⇡n n⇠alors
⇡ n alors
X n+1X⇠ ⇡?
n+1 ⇠ ⇡(y)
n+1 (y) K(x,
n+1= = K(x,
y)⇡ ny)⇡ n (x)dx
(x)dx
Une condition suffisante pour avoir une distribution stationnaire ⇡ est le principe “detailed balance”:
Detailed balance

Si ⇡ vérifie la condition:
K(x, y)⇡(x) = K(y, x)⇡(y) 8x, y

Alors (X n )n a une distribution stationnaire ⇡.

Interprétez cette équation en l’intégrant sur des régions A et B contenant x et y respectivement.

On obtient: P(Xn+1 2 A, Xn 2 B) = P(Xn+1 2 B, Xn 2 A)

“autant de particules vont de A à B que de B à A”

36
Markov Chains Monte-Carlo Dé nition
Algorithmes MCMC
def
Soit g une densité (non-normalisée). On veut simuler un échantillon x selon la densité ⇡ = Rg .
g

Soit (X n )n une chaı̂ne de Markov admettant une distribution stationnaire à densité ⇡.

Alors pour simuler x ⇠ ⇡ Il suffit de prendre x = X n avec n ! 1.

Pour générer un échantillon, il faut attendre la convergence de la chaîne

Et pour estimer une moyenne a posteriori, il faut plusieurs échantillons i.i.d !
Mais les itérés de la suite ne sont jamais i.i.d …
Théorème Ergodique

Soit (X n )n une chaı̂ne de Markov ergodique ayant pour distribution stationnaire ⇡. Alors:
n
X Z “Loi des grands
1 n!+1
nombres” pour des
f (X i ) ! E⇡ (f (X)) = f (x)⇡(x)dx
n i=1 échantillons non i.i.d

37
fi
Markov Chains Monte-Carlo Dé nition
Algorithmes MCMC
def
Soit g une densité (non-normalisée). On veut simuler un échantillon x selon la densité ⇡ = Rg .
g

Soit (X n )n une chaı̂ne de Markov admettant une distribution stationnaire à densité ⇡.

Alors pour simuler x ⇠ ⇡ Il suffit de prendre x = X n avec n ! 1.

Comment construire la suite (X n )n ?

1. Gibbs Sampling
Foundation of MCMC 2. Metropolis
3. Metropolis-Hastings (MH)
4. Langevin dynamics

Physics inspired algorithms 5. Hamiltonian Monte-Carlo (HMC)

6. No-U-Turn Sampler (NUTS)

Default algorithm in pyMC

38
fi
Markov Chains Monte-Carlo Gibbs

On veut souvent des échantillons suivant une loi a posteriori multivariée f (✓,µ,..., 2) . Peut-on gé
2 2
Où on sait générer ✓|µ, ..., et µ|✓, . . . , et ainsi de suite

2) . Peut-on générer un paramètre à la fois ?

Oui ! Grâce à l’algorithme de Gibbs.

39
Gibbs sampling En deux dimensions
Algorithme de Gibbs (2D)
On
On souhaite
souhaite simuler
simuler XX == (X
(X 1, ,
X X 2))22 RR2 2
d’une
d’une distribution
distribution jointe
jointe PP (X ,X ). . Soit
Soit xx 2 2
2 RR22
. .
On souhaite simuler X = (X 1
1 , X 2 ) 2 R 2
d’une distribution jointe
On souhaite simuler X = (X 1 , X 2 ) 2 R d’une distribution jointe P(X 1 ,X 2 ) . Soit x2 2 R .
2 2 P (X
(X 1
1
1,X
,X 2
2
2)
) . Soit x 2
2 2 R 2
2 .
1.
1. On
On sait
sait générer
générer xx 1 ⇠
⇠ X
X 1 avec
avec XX 2 fixée
fixée selon
selon la
la loi
loi conditionnelle
conditionnelle XX |X
1 |X 2 = = xx 2
1. On sait générer x 1 ⇠ X 1 avec X 2 fixée selon la loi conditionnelle
1. On sait générer x1 ⇠ X 1 avec X 2 fixée selon la loi conditionnelle X 1 |X 2 = x2 2
1 1 2 X 1
1 |X 2
2 = x 2

2.
2. On
On sait
sait générer
générer XX 2 avec
avec XX 1 fixée
fixée selon
selon lala loi
loi conditionnelle
conditionnelle XX 2 |X
|X 1 = = xx 1
2. On sait générer X 2 avec X 1 fixée selon la loi conditionnelle
2. On sait générer X 2 avec X 1 fixée selon la loi conditionnelle X 2 |X 1 = x1 1
2 1 X 2
2 |X 1
1 = x 1

Alors
Alors en
en partant
partant d’un
d’un quelconque
quelconque x x 2 22 R,
R, la
la suite
suite XX n définie
définie par
par ces
ces itérés
itérés est
est une
une chaı̂ne
chaı̂ne de
de
Alors en partant d’un quelconque x 2 2R,R, la suite X n définie par ces itérés est une chaı̂ne de
Alors
Markov
Markoven partant
avec
avec PP d’un
(X ,X quelconque
) comme
comme x 2
2
distribution
2
distribution la suite X
stationnaire.
stationnaire.n
ndéfinie par ces itérés est une chaı̂ne de
Markov avec P (X
(X 1
1,X
,X 2
2) comme distribution
Markov avec P(X 1 ,X 2 ) comme distribution stationnaire.
1 2 ) stationnaire.
Application (Ex 3. TD 1)
2 2 2
Soit Z|µ, ⇠ N (µ, ). µ et ont des densités 2 2
a priori indépendantes fµ et f 2 . Z|µ, ⇠ N (µ, )
1. Dessiner le graphe probabiliste du modèle.
2µ ⇠ f
2. Déterminez la loi a posteriori jointe (µ, 2
)|Z.
Z|µ, ⇠ Nµ (µ
2 2
2 2
Z|µ, ⇠ N (µ, )
3. Déterminez les lois conditionnelles µ| , Z et |µ, Z.
2
4. Quelles lois a priori fµ et f 2 devrait-on prendre pour Z|µ, ⇠ N⇠(µ,
f 2
2

avoir des lois conditionnelles usuelles ?

5. Implémentez l’algorithme de Gibbs pour simuler la loi a posteriori.

40
 Gibbs sampling En dimension d
Algorithme de Gibbs (général)

On souhaite simuler X = (X 1 , . . . , X d ) 2 Rd d’une distribution jointe P(X 1 ,...,X d ) . Soit x2 , . . . , xd 2 R.

1. Générer x1 ⇠ X 1 avec X 2 , . . . , X d fixés selon la loi conditionnelle X 1 |X 2 = x2 , . . . X d = xd

2. Générer (et mettre à jour) x2 ⇠ X 2 avec X 1 , X 3 , . . . , X d fixés selon la loi conditionnelle X 2 |X 1 = x1 , X 3 = x3 , . . . , Xd = xd
3. Générer (et mettre à jour) x3 ⇠ X 3 avec X 1 , X 2 , . . . , X d fixés selon la loi conditionnelle X 3 |X 1 = x1 , X 2 = x2 , . . . , Xd = xd
4. ...

La suite Xn définie par ces itérés est une chaı̂ne de Markov avec P(X 1 ,...,X n ) comme distribution stationnaire.

1. Simple et facile à mettre en oeuvre

Avantages:
2. E cace pour simuler les lois a posteriori d’un modèle hiérarchique

1. Nécessite de savoir simuler les lois conditionnelles

Inconvénients: 2. Convergence lente en grande dimension
3. Convergence très lente si les composantes sont corrélées
41
ffi
Intro to Metropolis

L’algorithme de Gibbs permet de découper le problème de simulation

d’une loi jointe en lois conditionnelles en 1D.

Et si la loi conditionnelle en 1D n’est pas usuelle ?

C’est là où intervient l’algorithme de Metropolis.

42
Metropolis Random Walk
2
Une chaîne de Markov simple est donnée par une marche aléatoire: X n+1 = X n + " avec " ⇠ N (0, )

Admet-elle une distribution stationnaire ?

Non ! Elle a une variance divergente (calcul simple)
Visuellement, c’est une suite qui “explore” l’espace:

43
Metropolis L’algorithme
On souhaite créer une chaîne de Markov avec la distribution stationnaire:

Quel est le problème avec cet algorithme ?

Idée de Metropolis:
Risque de rester coincé autour du maximum !
2
1. Explorer l’espace avec une marche aléatoire X n+1 = X n + " avec " ⇠ N (0, )

2. Utiliser la densité cible f pour accepter ou rejeter des “sauts”:

f (Xn+1 )
Accepter de sauter à X n+1 s’il est “meilleur” que X n : f (Xn+1 ) f (Xn ) avec probabilité f (Xn )

Sinon X n+1 = X n et refaire la marche aléatoire Comment l’implémenter en pratique ?

44
Metropolis L’algorithme
Algorithme de Metropolis (Gaussien)

Soit
Soitf ffune
Soit une densité
densitédede
unedensité probabilité.
deprobabilité. On
probabilité.On suppose
Onsuppose que
queXX
supposeque X
nn estdéjà
est déjàgénéré.
généré.XX
n est déjà généré. X
estdéfini
est définipar:
par:
n+1 est défini par:
n+1
n+1

2 2 2
1.1. Générery yy⇠⇠
Générer
[Link]́nérer ⇠NN (X
N(X
(Xn ,n
n
,
, 2)
) ) équivalent à y = Xn + " et générer " ⇠ N (0, )
2.2. Générer
[Link]́nérer u⇠⇠
Généreru u U([0,
⇠U([0, 1]).
U([0,1]).
1]). équivalent à accepter y avec probabilité min(1, f (y)
f (Xn ) )
f
f (y)
f (y)
(y)) )alors
3.3. Si
Si u u
< < min(1,
min(1,
3. Si u < min(1,f (X
f (X ) ) ) alors
X
alors X
X =
n+1
n+1
n+1
=
= y,
y,
y, sinon
sinon X
sinon X
X =
n+1
n+1
n+1
=
=XX
Xn .
n
n
.
.
f (Xnnn )

La suite (X n )n obtenue admet une distribution stationnaire donnée par la densité f .

45
Metropolis Metropolis in action

Exemple en deux dimensions:

Metropolis L’algorithme
Algorithme de Metropolis (Gaussien)

Soit
Soitf ffune
Soit une densité
densitédede
unedensité probabilité.
deprobabilité. On
probabilité.On suppose
Onsuppose que
queXX
supposeque X
nn estdéjà
est déjàgénéré.
généré.XX
n est déjà généré. X
estdéfini
est définipar:
par:
n+1 est défini par:
n+1
n+1

2 2
1.1. Générery yy⇠⇠
Générer
[Link]́nérer ⇠NN (X
N(X
(Xn ,n
n
,
, 2)
) )
2.2. Générer u⇠⇠
Généreru u
[Link]́nérer U([0,
⇠U([0, 1]).
U([0,1]).
1]).
f
f (y)
f (y)
(y)) )alors
3.3. Si
Si u u
< < min(1,
min(1,
3. Si u < min(1,f (X
f (X ) ) ) alors
X
alors X
X =
n+1
n+1
n+1
=
= y,
y,
y, sinon
sinon X
sinon X
X =
n+1
n+1
n+1
=
=XX
Xn .
n
n
.
.
f (Xnnn )

La suite (X n )n obtenue admet une distribution stationnaire donnée par la densité f .

Peut-on utiliser cet algorithme si on a uniquement accès à g / f ?

R f (y) g(y)
Oui ! f = Rg et la constante de normalisation g disparaı̂t dans le rapport =
g f (Xn ) g(Xn )

2 2
Quel est l’e↵et de ? contrôle le trade-o↵ “exploration-exploitation”

47
Metropolis L’algorithme
Algorithme de Metropolis (Gaussien)

Soit
Soitf ffune
Soit une densité
densitédede
unedensité probabilité.
deprobabilité. On
probabilité.On suppose
Onsuppose que
queXX
supposeque X
nn estdéjà
est déjàgénéré.
généré.XX
n est déjà généré. X
estdéfini
est définipar:
par:
n+1 est défini par:
n+1
n+1

2 2
1.1. Générery yy⇠⇠
Générer
[Link]́nérer ⇠NN (X
N(X
(Xn ,n
n
,
, 2)
) )
2.2. Générer u⇠⇠
Généreru u
[Link]́nérer U([0,
⇠U([0, 1]).
U([0,1]).
1]).
f
f (y)
f (y)
(y)) )alors
3.3. Si
Si u u
< < min(1,
min(1,
3. Si u < min(1,f (X
f (X ) ) ) alors
X
alors X
X =
n+1
n+1
n+1
=
= y,
y,
y, sinon
sinon X
sinon X
X =
n+1
n+1
n+1
=
=XX
Xn .
n
n
.
.
f (Xnnn )

La suite (X n )n obtenue admet une distribution stationnaire donnée par la densité f .

1. Il su t de savoir calculer une densité (même non-normalisée) f pour l’implémenter

Avantages:
2. Généralisable facilement en dimension > 1

1. La variance doit être assez petite pour éviter trop de rejets (random walk behavior)
Inconvénients: 2. La variance doit être assez grande pour bien explorer l’espace (distributions multimodales)
3. Le pourcentage d’acceptation tend vers 0 en grande dimension
ffi
Intro to Metropolis-Hastings

Qu’en est-il d’une densité cible discrète ?

Hastings généralise l’algorithme Metropolis où on peut choisir la

distribution qui explore l’espace (marche aléatoire).

49
Metropolis-Hastings L’algorithme
Hastings généralise l’algorithme de Metropolis en prenant une distribution de transition (exploration) quelconque:

Algorithme de Metropolis-Hastings

Soit f une densité de probabilité et q(y|x) une densité de transition (proposal) de x à y.

On suppose que Xn est déjà généré. X n+1 est défini par:

1. Générer y ⇠ q(.|Xn )
2. Générer u ⇠ U([0, 1]).
f (y)q(Xn |y)
3. Si u < min(1, f (Xn )q(y|Xn ) ) alors X n+1 = y, sinon X n+1 = Xn .

La suite (X n )n obtenue admet une distribution stationnaire donnée par la densité f .

Avantages: 1. Utile si la distribution cible f a un domaine borné (par ex. [0, 1], on peut choisir Q = Beta)
2. Utile si la distribution cible f est discrète (Q uniforme sur ensemble ni)
3. On peut utiliser une densité q non normalisée (Gaussienne tronquée)

Inconvénients: Identiques à ceux de l’algorithme de Metropolis

Remarque: si q est symétrique, on retrouve l’algorithme de Metropolis

fi
Monte-Carlo: récap des séances passées

def R
On est souvent confronté à calculer une quantité du type: I = '(x)f (x)dx.

Par exemple, avec:

Par exemple, avec:
Par exemple, avec:
[Link]
f = densité
exemple, non
avec: normalisée et ' = 1: I = constante de normalisation
1. f = densité non normalisée et ' = 1: I = constante de normalisation
1. f = densité non normalisée et ' = 1: I = constante de normalisation
2.
1. f
f =
= densité
densité normalisée
non et
normalisée' =
et 1
' =
A : I
1: = I Probabilité
= constante d’un
de évenement
normalisation
2. f = densité normalisée et ' = 1A : I = Probabilité d’un évenement
2. f = densité normalisée et ' = 1A : Ik = Probabilité d’un évenement
3.
2. ff =
= densité
densité normalisée
normalisée et
et '
' :=x17!: xIk = : IProbabilité
= moment d’und’ordre k de loi a posteriori.
évenement
3. f = densité normalisée et ' : x 7! A x : I = moment d’ordre k de loi a posteriori.
k
3. f = densité normalisée et ' : x 7! x : I = moment d’ordre k de loi a posteriori.
4. ff = densité ànormalisée
posteriorietnormalisée k ' : x 7! x: I = Estimateur de Bayes.
et
3.
4. f = densité ' : x !
7 x : I = moment d’ordre k de
= densité à posteriori normalisée et ' : x 7! x: I = Estimateur de Bayes. loi a posteriori.
4. f = densité à posteriori normalisée et ' : x 7! x: I = Estimateur de Bayes.
4. f = densité à posteriori normalisée et ' : x 7! x: I = Estimateur de Bayes.

51
Monte-Carlo: récap des séances passées
def R
On est souvent confronté à calculer une quantité du type: I = '(x)f (x)dx.
En pratique on utilise:

1. Un estimateur Monte-Carlo “brut” (Crude Monte-Carlo)

Pn Il su t de simuler (générer) les
Avec des échantillons X1 , . . . , Xn i.i.d ⇠ f on peut estimer I avec : Ib = 1
n i=1 '(Xi ). échantillons et calculer la
moyenne
Pour simuler les Xi , on utilise des transformations d’échantillons uniformes, Rejection sampling . . .

[Link]: Importance Sampling (préférentiel), Variables de contrôle … (à voir en Méthodes de simulation)

2. L’estimateur Markov-Chain Monte-Carlo

Pn
Avec une chaı̂ne de Markov ergodique (Xi )i à distribution stationnaire f , on peut estimer I avec : Ib = 1
n i=1 '(Xi ).

Pour simuler les Xi , on utilise:

Il faut véri er la convergence
1. Metropolis-Hastings (Exploration aléatoire + rejet) de la chaîne de Markov

2. Gibbs (loi jointe ! lois conditionnelles)

52
ffi
fi
MCMC diagnostics in 1D Visuals
Soit (X n )n 0 une chaı̂ne de Markov donnée par un algorithme MCMC.

b P
1
P
n n
L’estimateur Ib =I =
L’estimateur 1
n n i=1 '(X
i=1 '(X
i ) )
est
i est
de de bonne
bonne qualité
qualité si: si:

1. 1. Converge
Converge rapidement:
rapidement: lesles
XiXutilisés
i utilisés
ontont atteint
atteint le régime
le régime stationnaire.
stationnaire.
b P P n P P
petite variance: V(
b
2. L’estimateur a une petite variance: V(I) = = i=1i=1
2. L’estimateur a une I) n V('(X
V('(X i) 2
i) + + 2 i<ji<j Cov(X
Cov(X i ,jX
i, X ) j)
Régime stationnaire = X n ⇠ f :
la chaı̂ne semble “mélangée”,
pas de tendance,
Soit (X n )n 0 une chaı̂ne de Markov donnée par un algorithme MCMC. absence de “pattern”,
b 1
P n ne reste pas bloquée.
L’estimateur I = '(Xi ) est de bonne qualité si:
n i=1

1. Converge rapidement: les Xi utilisés ont atteint le régime stationnaire.

Soit (X ) n n 0 une chaı̂ne de Markov donnée par un a
hP P i
b = 2 n
2. L’estimateur a une petite variance: V(I) 1
n i=1 V('(Xi )) + 2 i<j Cov('(Xi ), '(Xj ))
Corr(X n , X n+k )

Faible variance = faibles corrélations Cov(X i , X j ).

Auto-corrélation avec lag k = Cov(X n , X n+k ).
Décroissante en fonction de k.
Empiriquement avec N échantillons:
AutoCorr(k) = Corr([X 0 , X 1 , . . . , X N k ], [X k , X 1+k , . . . , X N ])
Auto-corrélation avec lag k = Cov(X n , X n+k ). 53
MCMC diagnostics in 1D Visuals

“Trace plot”
Histogramme des échantillons

“Auto-correlation plot” Corr(X n , X n+k )

Auto-corrélation avec lag k = Cov(X n , X n+k ).

54
MCMC diagnostics in 1D Visuals

“Trace plot”
Histogramme des échantillons

En pratique, on “jette” la
première partie des
échantillons pour être dans le
régime stationnaire: période de
“chau e” / Burn-in.

“Auto-correlation plot”

lag k = Cov(X n , X n+k ).

La corrélation entre deux
échantillons consécutifs est ~ 0.3.
On peut prendre un échantillon sur
deux ou trois pour réduire les
corrélations: c’est le thinning

Auto-corrélation avec lag k = Cov(X n , X n+k ).

55
ff
MCMC diagnostics in 1D Visuals

Auto-corrélations excessives:

Corr(X n , X n+k )
Auto-corrélation avec lag k = Cov(X n , X n+k ).

Est-ce une bonne chaîne MCMC ? Raison principale ?

Comment expliquer ce comportement ? Metropolis: On n’explore pas assez l’espace

56
MCMC diagnostics in 1D Visuals
(autocorr inutile ici, on s’y
attendait vu le trace plot)

Corr(X n , X n+k )
Auto-corrélation avec lag k = Cov(X n , X n+k ).

Explications principales ?
Très probablement:
Est-ce une bonne chaîne MCMC ?
1. Erreur d’implémentation
Drift excessif !
2. Loi a priori extrêmement loin des données
3. Convergence très très lente
57
MCMC diagnostics in 1D Visuals
(autocorr inutile ici, on s’y
attendait vu le trace plot)

Corr(X n , X n+k )
Auto-corrélation avec lag k = Cov(X n , X n+k ).

Explications principales ?

Est-ce une bonne chaîne MCMC ? Très probablement:

1. Erreur d’implémentation
Aucun mélange visible, Chaîne bloquée !
2. Mauvais choix de la distribution de proposition
(exploration très limitée)

58
MCMC diagnostics in 1D Visuals

Corr(X n , X n+k )
Auto-corrélation avec lag k = Cov(X n , X n+k ).

En pratique:
1. Très rare d’avoir une loi a posteriori avec “un
Est-ce une bonne chaîne MCMC ? désert” entre les modes

Oui ! Mais la distribution cible est multimodale … 2. Si c’est le cas, il vaut mieux changer de modèle: Ici
par ex, ajouter une variable latente (si possible) Z =
0 ou 1 pour distinguer les modes

59
MCMC diagnostics in 1D Visuals
Lancer plusieurs chaînes avec des initialisations X0 di érentes:

Bon exemple Mauvais exemple

Comment dé nir “un

score” pour détecter
ces di érences ?

60
ff
fi
ff
MCMC diagnostics in 1D Metrics: R hat de Gelman-Rubin
Moyenne des moyennes X
Moyenne et m
n X
variance de (j) 2
Variance des moyennes “between” B = (X X)
chaque chaîne m 1 j
m
X
(1) (1) (1) 2 (1) 1 2 (j)
Chaı̂ne 1 : X1 , . . . , Xn !X , Moyenne des variances “within” W =
m j=1
..
.
(j) (j) 2 (j)
Chaı̂ne j : X1 , . . . , Xn(j) !X , Pourquoi n ? (j) 2
car V(X )= n
..
.
(m) (m) (m) 2 (m) Si les chaînes convergent B
Chaı̂ne m : X1 , . . . , X n !X , !0
vers la même moyenne: n
B ?
Sinon !b >0
n
“R hat” de Gelman-Rubin (1992) s
s 8
n 1 1
< = 1 si B ! 0 b def
n 1 2
n
1 2
within + n between
R =
b =def n W + nB !r n 2
R : within
W > 1 sinon
En pratique on veut R < 1.01
Détails en TD 61
MCMC diagnostics in 1D Metrics

s s
n 1 1 n 1 1
W+ B b = n W+ nB
= 1.025
def
b =
R = 1.005
def n n R
W W
Converging Not converging
62
MCMC diagnostics in 1D Metrics: E ective sample size (ESS)
On estime la moyenne d’une loi ⇡ avec deux estimateurs:

1. I0 : Moyenne de N0 échantillons i.i.d ⇠ ⇡

Quel est le meilleur estimateur ?
2. I1 : Moyenne de N1 échantillons ⇠ ⇡.

On suppose que V(I0 ) = V(I1 ). On peut en déduire que N0 N1 ou N0  N1 ?

ESS: E ective sample size (Hastings 1970)

On suppose
Étant donnée une MC X1 , . . . , XN1 , ESS correspond au nombre que V(I0 ) = V(I1 ). On peut en d
N0 tel que:
En général, sauf cas
N1
ESS = P1  N1 spéciaux (autocorr < 0)
1 + 2 t=1 AutoCorr(lag t) (variables antithétiques à voir
en Méthodes de simulation)
Preuve: il su t d’égaliser la variance des deux estimateurs en supposant le régime stationnaire atteint

Remarques

1. On interprète ESS comme le nombre d’échantillons “réellement” utilisés (en moins ~ 400)
2. On calcule une version améliorée qui prend en compte plusieurs chaînes (Gelman, 2013)

63
ff
ffi
ff
MCMC diagnostics in 1D Metrics
N = 5500

Corr(X n , X n+k )
ESS = 5149
échantillons “e caces”

Auto-corrélation avec lag k = Cov(X n , X n+k ).

N = 5500

Corr(X n , X n+k )
ESS = ESS = 145 << N
Notre estimation a la qualité
d’une moyenne i.i.d avec
145 échantillons seulement
Auto-corrélation avec lag k = Cov(X n , X n+k ).
Inconvénient: quanti e la qualité de l’estimation de la moyenne uniquement (ESS mean)

On peut dé nir d’autres ESS pour l’estimation des quantiles (ESS median, ESS Q95%..)

64
fi
ffi
fi
MCMC diagnostics in 1D Metrics
s
n 1 2 1 2
b def n within + n between
R = 2
within

Inconvénient:
On peut avoir R ~ 1 alors que les
chaînes sont très di érentes. Comment
construire un tel exemple ?

2
R⇡1, between ⇡0
Il su t d’avoir des chaînes di érentes
mais de mêmes moyennes !

s s
n 1 1 n 1 1
W + B W+ B
b =
R = 1.005
def n
ESS
n
= 4500 b =
R = 1.025
def n n
W W
Converging Not converging
65
ffi
ff
ff
MCMC diagnostics in 1D Metrics

“Rank-Normalization” (Vehtari et al, 2021)

(j) def (j)
1.
1. Trier
Trier toutes
toutes les
les observations
observations en
en calculant
calculant R
R
(j)
(j)
i =
def
def
= rank(X
rank(X
(j)
(j)
i )
) 2 [|1,
2 [|1, mn|]
mn|]
1. Trier toutes les observations en calculant Rii = rank(X ii ) 2 [|1, mn|]
(j)
(j)
(j) def X(j)
(j)
2.
2. Les
Les projeter
projeter sur
sur [0, 1]:
[0, 1]: u
u
(j)
(j) =
def
def
=
R
X
X
i
(j)
ii
i
2. Les projeter sur [0, 1]: uiii = mn
mn
mn
i

(j)
3. Si les chaı̂nes sont bien “mélangées” alors u (j)
(j) ⇠ U([0, 1])
3.
3. Si
Si les
les chaı̂nes
chaı̂nes sont
sont bien
bien “mélangées”
“mélangées” alorsalors uuii ⇠
i ⇠ U([0,
U([0, 1])
1])
def
(j) def 1 (j)(j)
4. Avec le théorème d’inversion Z (j)
(j) =
def 1 (u (j) ) ⇠ N (0, 1), :: Fct de répart. de N (0, 1).
4.
4. Avec
Avec lele théorème
théorème d’inversion
d’inversion Z Z ii =
i = 1
(uii )) ⇠
(u i ⇠NN (0,
(0, 1),
1), : Fct
Fct de
de répart.
répart. de
de N
N (0,
(0, 1).
1).
b (j) (j)
5.
5. Calculer
Calculer R
b
b classique sur Z (j)
(j) au lieu des X (j)
(j) .
5. Calculer R R classique
classique sur
sur Z
Z ii au
i au lieu
lieu des
des XX ii ..
i

Remarques:

1. Rhat plus sensible à la variabilité entre les chaînes

2. Moins sensible aux outliers
3. Cette transformation est utilisée pour calculer ESS également: “ESS bulk” (moyenne)
et “ESS tail” (quantile 0.95)
4. Méthode par défaut en Python (Arviz)
66
Récap de la séance passée
1.
1. Le
Le but
but des
des algorithmes
algorithmes MCMC
MCMC est
est de
de simuler
simuler des
des échantillons
échantillons suivant
suivant une
une loi
loi àà densité
densité connue
connue /
/ ff ..
1.
1. Le
Le but
but des
des algorithmes
algorithmes MCMC
MCMC est
est de
de simuler
simuler des
des échantillons
échantillons suivant
suivant une
une loi
loi àà densité
densité connue
connue /
/ ff ..
1.
1.
1.
2. Le
Le
Le
Un but
but
but des
des
des
algorithmealgorithmes
algorithmes
algorithmes
MCMC MCMC
MCMC
MCMC
crée une est
est
est de
de
chaı̂nede simuler
simuler
simuler
de Markov des
des
des échantillons
échantillons
échantillons
dont la suivant
suivant
suivant
distributions une
une
une loi
loi
loiàààdensité
asymptotique densité
densité connue
connue
estconnue
f /
// fff.. .
(stationnaire).
2.
1.
2.
1. Un
Le
Un
Le algorithme
but des
algorithme
but des MCMC
algorithmes
MCMC
algorithmes crée
MCMC
crée
MCMC une
une chaı̂ne
est
chaı̂ne
est de
de de
de Markov
simuler
Markov
simuler des
des dont
dont la
échantillons
la
échantillons distributions
suivant
distributions
suivant asymptotique
une loi
asymptotique
une loi à
à densité
densité est
est ff (stationnaire).
connue
connue / f
(stationnaire).
/ f .
.
2.
2. Un
Un algorithme
algorithme MCMC
MCMC crée
crée une
une chaı̂ne
chaı̂ne de
de Markov
Markov dont
dont la
la distributions
distributions asymptotique
asymptotique est
est f (stationnaire).
ff (stationnaire).
2.
2.
2.
3. Un
Un
Un
Ces algorithme
algorithme
algorithme
algorithmes MCMC
MCMC
MCMC
sont crée
crée
crée une
une
une
principalement chaı̂ne
chaı̂ne
chaı̂ne de
de
de
basés Markov
Markov
Markov
sur une dont
dont
dont la
la
explorationla distributions
distributions
distributions
aléatoire asymptotique
asymptotique
asymptotique
et un mécanisme est
est
est ff
de (stationnaire).
(stationnaire).
(stationnaire).
rejet.
3.
2.
3.
2. Ces
Un
Ces
Un algorithmes
algorithme
algorithmes
algorithme sont
MCMC
sont
MCMC principalement
crée une
principalement
crée une chaı̂ne
chaı̂nebasés
de
basés
de sur
Markov
sur
Markovune
une exploration
dont
exploration
dont la
la aléatoire
distributions
aléatoire
distributions et
et un
un mécanisme
asymptotique
mécanisme
asymptotique est
est de
f
de
f rejet.
(stationnaire).
rejet.
(stationnaire).
3.
3. Ces
Ces algorithmes
algorithmes sont
sont principalement
principalement basés
basés sur
sur une
une exploration
exploration aléatoire
aléatoire et
et un
un mécanisme
mécanisme de
de rejet.
rejet.
3.
3.
3.
4. Ces
Ces
Ces
En algorithmes
algorithmes
algorithmes
pratique, on sont
sont
sont
jette principalement
principalement
principalement
les k premiers basés
basés
basés sur
sur
sur
échantillons une
une
une
pourexploration
exploration
exploration
être dans aléatoire
aléatoire
aléatoire
le régime et
etet un
un
un mécanisme
mécanisme
mécanisme
stationnaire: de
de
de
burn-in rejet.
rejet.
rejet.
/ tuning.
4.
3.
4.
3. En
Ces
En
Ces pratique,
pratique, on
algorithmes
on
algorithmes jette
sont
jette
sont les
les k
k premiers
principalement
premiers
principalement échantillons
basés sur
échantillons
basés sur pour
une
pour
une être
exploration
être
exploration dans
dans le
le régime
aléatoire
régime
aléatoire stationnaire:
et un mécanisme
stationnaire:
et un mécanisme burn-in
de
burn-in
de /
/ tuning.
rejet.
tuning.
rejet.
4.
4. En
En pratique,
pratique, on
on jette
jette les
les k
k premiers
premiers échantillons
échantillons pour
pour être
être dans
dans le
le régime
régime stationnaire:
stationnaire: burn-in
burn-in /
/ tuning.
tuning.
4.
4.
4.
5. En
En
En
En pratique,
pratique,
pratique,
pratique, on
on
laon jette
jette
jette
chaı̂ne les
les
lesk
kk premiers
premiers
premiers
obtenue n’est échantillons
échantillons
échantillons
pas forcément pour
pour
pour être
être
être
convergente, dans
dans
dans le
lelerégime
régime
pourrégime
cela stationnaire:
stationnaire:
stationnaire:
on doit burn-in
burn-in
e↵ectuerburn-in
un /
// tuning.
tuning.
tuning.
diagnostic de convergence.
5.
4.
5.
4. En
En
En
En pratique,
pratique,
pratique,
pratique, la
laon
on chaı̂ne
jette
chaı̂ne
jette obtenue
les k
obtenue
les k n’est
premiers
n’est
premiers pas
pas forcément
échantillons
forcément
échantillons pour
pourconvergente,
être
convergente,
être dans
dans pour
le
pour
le cela
régime
cela
régime on
on doit
stationnaire:
doit
stationnaire: e↵ectuer
e↵ectuer un
burn-in
un
burn-in diagnostic
/ tuning.
diagnostic
/ tuning. de
de convergence.
convergence.
5.
5. En
En pratique,
pratique, la
la chaı̂ne
chaı̂ne obtenue
obtenue n’est
n’est pas
pas forcément
forcément convergente,
convergente, pour
pour cela
cela on
on doit
doit e↵ectuer
e↵ectuer un
un diagnostic
diagnostic de
de convergence.
convergence.
5.
5.
5.
6. En
En
En
Chaquepratique,
pratique,
pratique, la
la
chaı̂ne la chaı̂ne
chaı̂ne
chaı̂ne
dépend obtenue
obtenue
obtenue
de son n’est
n’est
n’est pas
pas
pas
initialisation, forcément
forcément
forcément
pour celaconvergente,
convergente,
convergente,
on lance pour
pour
pour
plusieurs cela
cela
cela on
on
on
chaı̂nesdoit
doit
doit e↵ectuer
e↵ectuer
pour e↵ectuer
comparer un
unun diagnostic
diagnostic
diagnostic
leur de
de
de convergence.
convergence.
convergence.
distribution asymptotique.
5.
6.
5.
6. En
Chaque
En
Chaque pratique,
chaı̂ne
pratique,
chaı̂ne la
la chaı̂ne
dépend
chaı̂ne
dépend obtenue
de son
obtenue
de son n’est
n’est pas
initialisation,
pas
initialisation, forcément
pour
forcément
pour cela
cela convergente,
on lance
convergente,
on lance pour
plusieurs
pour
plusieurs cela
cela on
chaı̂nes
on
chaı̂nes doit
doitpour
pour e↵ectuer
comparer
e↵ectuer
comparer un
un diagnostic
leur
diagnostic
leur de
de convergence.
distribution asymptotique.
convergence.
distribution asymptotique.
6.
6. Chaque
Chaque chaı̂ne
chaı̂ne dépend
dépend de
de son
son initialisation,
initialisation, pour
pour cela
cela on
on lance
lance plusieurs
plusieurs chaı̂nes
chaı̂nes pour
pour comparer
comparer leur
leur distribution
distribution asymptotique.
asymptotique.
6.
6.
7.
7. Chaque
Chaque
Visuellement,
Visuellement,chaı̂ne
chaı̂neles
lesdépend
dépend
chaı̂nes
chaı̂nesde
de son
son
doivent
doivent initialisation,
initialisation,
être
être bien
bien pour
pour
mélangées
mélangéescela
cela on
et
et on lance
lance
des
des plusieurs
plusieurs
autocorrélations
autocorrélations chaı̂nes
chaı̂nes
faibles
faibles pour
pour (<
(<comparer
comparer
0.1).
0.1). leur
leur distribution
distribution asymptotique.
asymptotique.
6.
7.
6.
7. Chaque
Visuellement,
Chaque
Visuellement, chaı̂ne
chaı̂neles
les dépend
chaı̂nes
dépend
chaı̂nes de
de son
doivent
son
doivent initialisation,
être bien
initialisation,
être bien pour
mélangées
pour
mélangées cela
cela et
et on
on lance
des
lance
des plusieurs
autocorrélations
plusieurs
autocorrélations chaı̂nes
faibles
chaı̂nes
faibles pour
pour (<
(< comparer
0.1).
comparer
0.1). leur
leur distribution
distribution asymptotique.
asymptotique.
7.
7. Visuellement,
Visuellement, les
les chaı̂nes
chaı̂nes doivent
doivent être
être bien
bien mélangées
mélangées et
et des
des autocorrélations
autocorrélations faibles
faibles (<
(< 0.1).
0.1).
7.
7.
8. Visuellement,
Visuellement,
On calcule le les
R̂ les
dechaı̂nes
chaı̂nes doivent
doivent
Gelman-Rubin être
être
qui bien
bien
mesuremélangées
mélangées
la etet
variabilité des
des autocorrélations
autocorrélations
entre ces chaı̂nes. faibles
faibles
Il faut (<
(< 0.1).
avoir0.1).
un R̂ < 1.01.
7.
8.
7.
8. Visuellement,
On calcule
Visuellement,
On calcule le
le R̂
R̂ les
de
les
de chaı̂nes
chaı̂nes doivent
Gelman-Rubin
doivent
Gelman-Rubin être
qui
être
qui bien
mesure
bien
mesure mélangées
la
mélangées
la et
variabilité
et
variabilité des
des autocorrélations
entre ces chaı̂nes.
autocorrélations
entre ces chaı̂nes. faibles
Il faut
faibles
Il faut (<
avoir
(<
avoir0.1).
un
0.1).
un R̂
R̂ <
< 1.01.
1.01.
8.
8. On
On calcule
calcule le
le R̂
R̂ de
de Gelman-Rubin
Gelman-Rubin qui
qui mesure
mesure la
la variabilité
variabilité entre
entre ces
ces chaı̂nes.
chaı̂nes. Il
Il faut
faut avoir
avoir un
un R̂
R̂ <
< 1.01.
1.01.
8.
8.
9.
9. On
On
On
On calcule
calcule
calcule
calcule le
le R̂
l’ESS
l’ESS R̂ dede Gelman-Rubin
Gelman-Rubin
(e↵ective
(e↵ective sample
sample qui
qui
size)
size) mesure
mesure
qui
qui la
donne
donne la variabilité
variabilité
une
une mesure
mesure entre
entre
du
du ces
ces chaı̂nes.
chaı̂nes.
nombre
nombre Il Il faut
faut
d’échantillons
d’échantillons avoir
avoir un
un R̂
R̂
“réllement
“réllement << 1.01.
1.01.
utilisés”.
utilisés”.
8.
9.
8.
9. On
On
On
On calcule
calcule
calcule
calcule le
l’ESS
le
l’ESS R̂
R̂ de
de Gelman-Rubin
(e↵ective sample
Gelman-Rubin
(e↵ective sample qui
size)
qui
size) mesure
qui
mesure
qui donne
donne la
la variabilité
une mesure
variabilité
une mesure entre
du
entre
du ces
ces chaı̂nes.
nombre
chaı̂nes.
nombre Il
Il faut
d’échantillons
faut
d’échantillons avoir
avoir un
un R̂
“réllement
R̂
“réllement <
< 1.01.
utilisés”.
1.01.
utilisés”.
9.
9. On
On calcule
calcule l’ESS
l’ESS (e↵ective
(e↵ective sample
sample size)
size) qui
qui donne
donne une
une mesure
mesure du
du nombre
nombre d’échantillons
d’échantillons “réllement
“réllement utilisés”.
utilisés”.
9.
9.
10. On
On
ESS calcule
calcule
bulk l’ESS
l’ESS
donne (e↵ective
le (e↵ective
nombre sample
sample size)
d’échantillonssize) qui
qui donne
pourdonne une
une
l’estimateur mesure
mesurede du
du
la nombre
nombre
moyenne. d’échantillons
d’échantillons “réllement
“réllement utilisés”.
utilisés”.
9.
10.
9.
10. On
ESS
On
ESS calcule
bulk
calcule
bulk l’ESS
donne
l’ESS
donne le
le (e↵ective
nombre
(e↵ective
nombre sample
d’échantillons
sample
d’échantillonssize)
size) qui
pour
qui
pourdonne
donne une
l’estimateur
une
l’estimateur mesure
mesurede
de du
la
du
la nombre
moyenne.
nombre
moyenne. d’échantillons
d’échantillons “réllement
“réllement utilisés”.
utilisés”.
10.
10. ESS
ESS bulk
bulk donne
donne le
le nombre
nombre d’échantillons
d’échantillons pour
pour l’estimateur
l’estimateur de
de la
la moyenne.
moyenne.
10.
10.
11.
11. ESS
ESS
ESS
ESS bulk
bulk
tail
tail donne
donne
donne
donne le
lele
le nombre
nombre
nombre
nombre d’échantillons
d’échantillons
d’échantillons
d’échantillons pour
pour
pour
pour l’estimateur
l’estimateur
l’estimateur
l’estimateur de
de
des
des la
la moyenne.
moyenne.
quantiles
quantiles 0.05
0.05 et
et 0.95
0.95 (le
(le min
min des
des deux
deux ESS).
ESS).
10.
11.
10.
11. ESS
ESS
ESS
ESS bulk
tail
bulk
tail donne
donne
donne
donne le
le le
le nombre
nombre
nombre
nombre d’échantillons
d’échantillons
d’échantillons
d’échantillons pour
pour
pour
pour l’estimateur
l’estimateur
l’estimateur
l’estimateur de
des
de
des la
la moyenne.
quantiles
moyenne.
quantiles 0.05
0.05 et
et 0.95
0.95 (le
(le min
min des
des deux
deux ESS).
ESS).
11.
11. ESS
ESS tail
tail donne
donne le
le nombre
nombre d’échantillons
d’échantillons pour
pour l’estimateur
l’estimateur des
des quantiles
quantiles 0.05
0.05 et
et 0.95
0.95 (le
(le min
min des
des deux
deux ESS).
ESS).
11.
11.
12. ESS
ESS
En tail
tail
pratiquedonne
donne
il le
fautlenombre
nombre
des ESS d’échantillons
d’échantillons
> 400 pour pour
pour
avoir l’estimateur
l’estimateur
une estimation des
desquantiles
de quantiles
bonne 0.05
0.05 et
qualité. et 0.95
0.95 (le
(le min
min des
des deux
deux ESS).
ESS).
11.
12.
11.
12. ESS
En
ESS
En tail
pratique
tail
pratique donne
il
donne
il faut
fautle
le nombre
des ESS
nombre
des ESS d’échantillons
> 400 pour
d’échantillons
> 400 pour pour
avoir
pour
avoir l’estimateur
une estimation
l’estimateur
une estimation des
de
des
de quantiles
bonne
quantiles
bonne 0.05
qualité.
0.05
qualité. et
et 0.95
0.95 (le
(le min
min des
des deux
deux ESS).
ESS).
12.
12. En
En pratique
pratique il
il faut
faut des
des ESS
ESS >
> 400
400 pour
pour avoir
avoir une
une estimation
estimation de
de bonne
bonne qualité.
qualité.
12.
12.
13.
13. En
En
Si
Si le
lepratique
pratique il
diagnostic
diagnosticil faut
faut
est
est des
des ESS
ESS
négatif,
négatif, >
il
il > 400
400
faut
faut pour
pour
identifier
identifieravoir
avoir
la
la une
une
cause
cause estimation
estimation
et
et changer
changer de
de bonne
les
les bonne qualité.
qualité.
paramètres
paramètres ou
ou le
le sampler.
sampler.
12.
13.
12.
13. En
Si
En
Si le
le pratique
diagnostic
pratique
diagnosticil
il faut
est
faut
est des
des ESS
négatif,
ESS
négatif, il
il >
> 400
faut
400
faut pour
identifier
pour
identifier avoir
la
avoir
la une
cause
une
cause estimation
et
estimation
et changer
changer de
deles
les bonne
bonne qualité.
paramètres
qualité.
paramètres ou
ou le
le sampler.
sampler.
13.
13. Si
Si le diagnostic est négatif, il faut identifier la cause et changer les paramètres ou le sampler.
13.
13. SiSile
lelediagnostic
diagnostic
diagnostic est
est
est négatif,
négatif,
négatif, il
ililfaut
faut
faut identifier
identifier
identifier la
lalacause
cause
cause et
etetchanger
changer
changer les
les
les paramètres
paramètres
paramètres ou
ouou le
lelesampler.
sampler.
sampler.
13.
13. Si le diagnostic est négatif, il faut identifier la cause et changer les paramètres ou le sampler.
Si le diagnostic est négatif, il faut identifier la cause et changer les paramètres ou le sampler.

67
MCMC diagnostics in 1D MCSE
Monte-Carlo Standard Error (MCSE): L’écart type d’un estimateur Monte-Carlo
def
On suppose le régime stationnaire atteint Xi ⇠ f . On pose = V(Xi )
n 2 3
1 X n
X X 2
L’estimateur de la moyenne est: Xi Sa variance: 1 4 5 =
n V(X i ) + 2 Cov(X i , X j )
i=1 n2 i=1 i<j
n ESS
" n 1
#
1 X X
ˆ = 2 +2 1 Cov(X
Xn i , Xi+k )
n ˆi=1 2
MCSE = p = k=1(X X̄)
i
nESS n 1 i=1

MCSE de la moyenne

Avec le même principe, on peut dé nir le MCSE de la médiane, de l’écart type lui même etc …

68
fi
MCMC diagnostics in 1D HDI: High density intervals
High density interval (HDI) ou Credible Interval (CI):
équivalent de l’intervalle de con ance en statistiques bayésiennes

Comment dé nir un intervalle à haute densité de niveau 90% ?

Un intervalle I tel que Pf (X 2 I) 0.9

R
Un intervalle I tel que I
f 0.9

On peut en trouver plusieurs !

Lequel faut-il choisir ?

HDI de niveau 1 ↵ = Le plus petit intervalle I telle que Pf (X 2 I) 1 ↵

Comment peut-on les comparer aux intervalles de con ance d’un point de vue statistique ?

69
fi
fi
fi
Application (Python)
Application (MCMC Diagnostics)
Soit Z|µ, 2 ⇠ N (µ, 2 ). µ et 2 ont des densités
2 2
a priori indépendantes fµ et f 2 . Z|µ, ⇠ N (µ, )
1. Dessiner le graphe probabiliste du modèle.
2µ ⇠ f
2. Déterminez la loi a posteriori jointe (µ, 2
)|Z. Z|µ, ⇠ Nµ (µ,
2 2
3. Déterminez les lois conditionnelles µ| 2
, Z et 2
|µ, Z. Z|µ, ⇠ N (µ, )
4. Quelles lois a priori fµ et f 2 devrait-on prendre pour 2 2
avoir des lois conditionnelles usuelles ?
Z|µ, ⇠ N⇠(µ,
2
f 2

5. Implémentez l’algorithme de Gibbs pour simuler la loi a posteriori.

Simuler 4 chaînes en numpy (Gibbs) et faîtes le diagnostique avec Arviz:

70
II Méthodes de Monte-Carlo
1. Introduction
2. Markov Chain Monte-Carlo (MCMC)
3. Algorithmes MCMC avancés

71
MCMC: algorithmes avancés Langevin Dynamics
Des particules d’un uide dont le mouvement est expliqué par:
d 2 Xt dXt p
1. Des forces déterministes
m 2 (gravitation)
= + F (Xt ) + 2 kT ⇠(t)
dt dt
2. Forces aléatoires (chaleur thermique) G(Xt )
d 2 Xt dXt p
3. Frottements / amortissement
m 2 = + F (Xt ) + 2 kT ⇠(t)
dt dt
Loi de Newton généralisée:
The Langevin Equation (1908)
2
d Xt dXt
m 2 = + F (Xt ) + G(Xt )
dt dt
p
2
d Xt
m 2 ⇡0 F (Xt ) = rU (Xt ) G(Xt ) = 2 kT ⇠t ⇠t ⇠ N (0, 1)
dt
Dé nition de l’énergie potentielle U Moyenne des forces aléatoires modélisée
La masse d’une particule
par une Gaussienne (thm. central)
est négligeable
s k: constante de Boltzmann

Ainsi:
dXt 1 2kT
= rU (Xt ) + ⇠t
dt
72
fi
fl
MCMC: algorithmes avancés Langevin Dynamics
s
dXt 1 2kT
= rU (Xt ) + ⇠t
dt
✓ ◆
U (x)
À l’équilibre, la distribution de ces particules est celle de Boltzmann: f (x) / exp
kT

U (x) rU (x)
) log (f (x)) = + cte ) r log (f (x)) =
kT kT
s
dXt kT 2kT
= rlog (f (Xt )) + ⇠t
dt Coe cient de di usion
d’Einstein
Stochastic di erential dXt p def kT
= Drlog (f (Xt )) + 2D⇠t D =
equation (SDE) dt

Pouvez-vous deviner un moyen de générer des échantillons suivant une densité f ?

Si on arrive à résoudre cette équation différentielle alors on aura un processus ~ f

73
ffi
ff
ff
MCMC: algorithmes avancés Langevin Dynamics
p Comment résoudre cette SDE
dXt = Drlog (f (Xt )) dt + 2D⇠t dt numériquement ?
Z t+h Z t+h p Z t+h
⇠t : processus aléatoire ⇠ N
dXs = D rlog (f (Xs )) ds + 2D ⇠s ds ⇠t : processus aléatoire ⇠ N (0, 1)
t t t ⇠t : processus aléatoire ⇠ N (0, 1)
p
⇡ “Somme de h variables p
N (0, 1) i.i.d”
Avec h ! 0: ⇡ “Somme de h variables N (0,=1) i.i.d
hN (0, 1) = h⇠t
en t”

p
p
p
XX
Xt+h
t+h
t+h
t+h XX
t tt=
X =Dhrlog
t= Dhrlog(f
Dhrlog (f(X
(f (Xttt))
(X )) dt
dt + 2Dh⇠tttt
+ 2Dh⇠
2Dh⇠

Avec D = 1 et " = h ⇡ 0, on définit la suite: C’est la méthode d’Euler

p
Xn+1 = Xn + "rlog (f (Xn )) + 2"⇠n ⇠n ⇠ N (0, 1)

La chaı̂ne de Markov (Xn )n a pour distribution stationnaire une densité / f .

74
MCMC: algorithmes avancés Langevin Dynamics
Unadjusted Langevin Algorithm (ULA)
p
Xn+1 = Xn + "rlog (f (Xn )) + 2"⇠n ⇠n ⇠ N (0, 1)
Interpréter cette suite. Quel est le lien avec Metropolis ?
1. Dirige la suite vers les régions à haute densité en explorant l’espace
2. Si on enlève le terme du gradient, on retrouve l’étape d’exploration de Metropolis
Pas d’étape de rejet: risque de non-convergence due aux erreurs de discrétisation

Metropolis Adjusted Langevin Algorithm (MALA)

1
Soit f une densité de probabilité et q(y|x) / exp( 4" ky x "r log f (x)k2 ). On suppose que Xn
est déja généré. X n+1 est défini par:
1. Générer y = Xn + "r log f (Xn ) + N (0, 2").

2. Générer u ⇠ U([0, 1]).

f (y)q(Xn |y)
3. Si u < min(1, f (Xn )q(y|Xn ) ) alors X n+1 = y, sinon X n+1 = Xn .

La chaı̂ne (X n )n obtenue admet une distribution stationnaire donnée par la densité f .

75
MCMC: algorithmes avancés Metropolis vs MALA

(animated)
76
MCMC: algorithmes avancés Metropolis vs MALA

Quand faut-il utiliser MALA au lieu de Metropolis ?

77 (animated)
MCMC: algorithmes avancés Metropolis vs MALA

La densité de Boltzmann donne une intuition très importante

✓ ◆
U (x)
f (x) / exp Comment interpréter “l’énergie potentielle” U ?
kT
haute énergie potentielle U(x) = hauteur (gravité)

Densité f Énergie potentielle -log(f)

78
MCMC: algorithmes avancés Hamiltonian Dynamics
En grande dimension, les algorithmes Metropolis / MALA explorent ine cacement l’espace.
Idée: au lieu d’explorer avec une position aléatoire x, simuler une vitesse initiale aléatoire.

En physique, on parle du moment: p~ = m~v

Hamiltonian Monte-Carlo (Duane 1987, Neal 1996)

1. Simuler un vecteur de moment p’ (direction et vitesse)

2. Déterminer la trajectoire
3. Suivre la trajectoire et s’arrêter (nouveau x’)
4. Accepter ou rejeter (x’, p’)
5. Répéter

(animated)
79

ffi
MCMC: algorithmes avancés Hamiltonian Dynamics
Hamiltonian Monte-Carlo (Duane 1987, Neal 1996)

1. Simuler un vecteur de moment p’ (direction et vitesse)

2. Déterminer la trajectoire
3. Suivre la trajectoire et s’arrêter (nouveau x’)
4. Accepter ou rejeter (x’, p’)
5. Répéter

Avec un moment aléatoire p', comment obtenir la trajectoire ?

Principe de conservation de l’Énergie totale = Énergie potentielle + Énergie cinétique

Énergie totale = l’opérateur Hamiltonien: H(x, p) = U (x) + K(p)

dH
Conservation de l’énergie dans le temps: (x, p) = 0
dt
def
Pour trouver la trajectoire, on discrétise l’équation di érentielle à p’ xé avec un pas " = t
On s’arrête après un nombre de pas égal à L En pratique, L et " sont difficiles à choisir
80
ff
fi
MCMC: algorithmes avancés Hamiltonian Dynamics

No-U-Turn-Sampler (NUTS) (Ho man and Gelman, 2014)

HMC
HMC avec:
avec:
HMC
HMC avec:
avec:
1.
1. La
La variance
variance d’exploration
d’exploration (pour
(pour générer
générer les
les moments)
moments) est
est adaptée
adaptée pendant
pendant le
le burn-in.
burn-in.
1.
1. La variance d’exploration (pour générer les moments) est adaptée pendant le burn-in.
La variance d’exploration (pour générer les moments) est adaptée pendant le burn-in.
2.
2. LL est
est optimisé
optimisé pour
pour éviter
éviter un
un U-turn
U-turn (demi-tour).
(demi-tour).
2.
2. L est optimisé pour éviter un U-turn (demi-tour).
L est optimisé pour éviter un U-turn (demi-tour).
3.
3. "" est
est modifié
modifié de
de façon
façon adaptative
adaptative pour
pour un
un taux
taux d’acceptation
d’acceptation cible.
cible.
3.
3. " est modifié de façon adaptative pour un taux d’acceptation cible.
" est modifié de façon adaptative pour un taux d’acceptation cible.

1. NUTS est l’algorithme “état-de-l’art” pour simuler en grande dimension (100-1000)

2. Aucun paramètre à choisir
3. Nécessite le gradient du log de la densité (distributions continues uniquement)
4. Par défaut dans PyMC pour les variables continues

81
ff
MCMC: algorithmes avancés Hamiltonian Dynamics

Barrière d’énergie

Est-il possible qu’une chaîne simulée par NUTS ne découvre qu’un seul mode de cette distribution ?

Oui ! Pour sortir de l’un des deux puits, il faut simuler un moment (vitesse) très grand

82
Algorithmes MCMC Combining MCMC algorithms

On veut simuler une loi a posteriori (↵, m|data) avec ↵ continue et m discrète. Comment faire ?

Algorithme de Gibbs:

Simuler ↵n ⇠ ↵|mn avec une itération de NUTS / MALA

Simuler mn+1 ⇠ m|↵n avec une itération de Metropolis-Hastings

Géré automatiquement par PyMC: à chaque variable (ou groupe de variables) est associé le meilleur
algorithme (exemple en TP).

83
Algorithmes MCMC Résumé

Algorithme Méthode Utilité

Lois conditionnelles simples

Gibbs Alterne entre les lois conditionnelles
Alterner entre des variables discrètes / continues

Marche aléatoire avec densité de proposition

Metropolis Densité non di érentiable symétrique
symétrique

Marche aléatoire avec densité de proposition Densité non di érentiable (discrète) et asymétrique
Metropolis-Hastings
quelconque comme la LogNormal, Gamma

Metropolis Adjusted
Marche aléatoire Gaussienne guidée par le
Langevin Algorithm Variable continue en faible dimension
gradient de la densité
(MALA)

Hamiltonian Monte- Exploration avec des moments aléatoires et des

Variable continue en grande dimension
Carlo (HMC) trajectoires à énergie totale constante

No-U-Turn Sampler Variable continue en grande dimension, pas de

HMC avec des trajectoires optimisées
(NUTS) paramètres à xer manuellement

84
fi
ff
ff
Chapitre 3. Applications et thématiques avancées

1. Modèles Bayésiens hiérarchiques (Assurance / Biostats)

2. Bayesian Machine learning

85
Bibliography

86