Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Chapitre 5 : Monté Carlo par Chaîne de
Markov - Monte Carlo Markov Chain
Dr. Nawel Arrar
3 ème année ENSIA
Janvier 2024
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 1/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Introduction
Il s’agit d’une méthode pour générer une suite fX1 , X2 , ...g de
variables aléatoires pouvant être interprétée comme constituant
une trajectoire d’une chaîne de Markov sur M états, possédant la
distribution stationnaire π = (π 1 , , π M ) avec 8i, π i 0 et
∑i π i = 1.
Cette méthode s’avère particulièrement riche en applications qui
vont de la génération d’états de systèmes physiques en équilibre
thermodynamique à l’obtention de points aléatoires dans des
régions compliquées.
Elle peut-être utilisée dans l’intégration par la méthode
Monté-Carlo lorsque la convergence de l’estimateur bIn vers I est
valable si les variables Xn ne sont plus indépendantes, mais
constituent par exemple une chaîne de Markov de loi stationnaire
f . On parle alors de Monte-Carlo par chaîne de Markov (MCMC),
dont l’algorithme de Metropolis-Hastings est l’exemple le plus
connu.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 2/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Introduction
Cette méthode, appelée méthode de Metropolis, a également
mené au développement d’heuristiques puissantes d’optimisation
combinatoire et continue. En statistique, cette approche, connue
sous l’acronyme MCMC (Markov Chain Monte Carlo), est
également très fréquemment employée.
Basée sur une analogie avec le procédé du recuit (courant en
métallurgie, il consiste à réchau¤er, puis refroidir très lentement un
alliage, pour lui permettre d’atteindre un état minimisant son
énergie interne), la méthode de Cerny-Kirkpatrick est connue sous
le nom de recuit simulé. Elle est également très utilisée dans la
phase d’apprentissage dans certains réseaux de neurones arti…ciels,
tel que cela a été traité en particulier dans le livre de Aarts et
Korst sur les machines de Boltzmann.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 3/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Propriétés des chaînes de Markov
De…nition
Soit π = (π 1 , , π M ) une distribution de probabilité et (Xn )
une chaîne de Markov de matrice de transition P. On dira que π
est réversible pour la chaîne si et seulement si elle véri…e les
équations d’équilibre détaillé :
π i pij = π j pji , 81 i, j M. (1)
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 4/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Algorithme de Metropolis
Nous présentons en premier lieu l’algorithme proposé par
Metropolis et ses co-auteurs en 1953 avant de passer à sa
généralisation par Hastings en 1970. On veut simuler une chaîne
de Markov de loi stationnaire π, celle-ci n’étant éventuellement
connue qu’à une constante de normalisation près :
π uj
8j 2 f1, , Mg π j = C π uj = , (2)
∑M u
i =1 π i
où π u signi…e “π unnormalized”. Ce problème est tout sauf
arti…ciel : il est par exemple récurrent en statistique bayésienne et
en physique statistique. Nous supposerons pour simpli…er, que
π j > 0 pour tout j.
Pour construire cette chaîne, il nous faut spéci…er une matrice de
transition P ayant π comme loi stationnaire. A l’instar des
méthodes de rejet et d’échantillonnage préférentiel.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 5/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Algorithme de Metropolis
Soit donc Q = [qij ] 1 i, j M une matrice de transition telle
que :
- symétrie : pour tout couple (i, j ), on a qij = qji ;
- partant de tout état i, on sait simuler facilement
selon la loi qi = [qi1 , ..., qiM ].
Algorithme
1 partant de Xn = i, simuler Y = j selon la loi
qi . = [qi1 , ..., qiM ] ;
πj
2 calculer le rapport d’acceptation rij = πi ;
3 tirer une loi uniforme U U [0, 1], et poser
Y = j si U rij ,
Xn +1 =
Xn = i sinon.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 6/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Algorithme de Metropolis
Proposition (Réversibilité de Metropolis) Si la matrice de
transition Q est irréductible, alors la chaîne (Xn ) produite par
l’algorithme de Metropolis est irréductible et π est réversible pour
cette chaîne. En particulier, π est son unique distribution
d’équilibre.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 7/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
Dans l’algorithme précédent, la dynamique dé…nie par la matrice
de transition Q était supposée agir symétriquement sur l’espace
d’états, à savoir que
qij = qji 8(i, j ) 2 E E.
Dans les applications, cette condition n’est pas nécessaire. Pour
s’en émanciper, Hastings a proposé une démarche permettant de
partir d’une matrice Q d’une chaîne de Markov irréductible. Mais,
il su¢ t de tenir compte de l’éventuelle dissymétrie de Q dans le
rapport rij . C’est ainsi que Hastings a généralié de l’algorithme de
Metropolis.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 8/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
Algorithme de Metropolis-Hastings
Pour cette généralisation étant simplement de supposer que :
- partant de tout état i, on sait simuler facilement
selon la loi qi. = [qi1 , ..., qiM ].
- pour tout couple (i, j ), on a qij > 0 implique qji > 0 ;
- pour tout couple (i, j ) tel que qij > 0, on sait
πq
calculer πji qijji .
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 9/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
Algorithme de Metropolis-Hastings
Algorithme
1 partant de Xn = i, simuler Y = j selon la loi
qi. = [qi1 , ..., qiM ] ;
2 calculer le rapport de Metropolis-Hastings
π j qji
rij = ,
π i qij
3 tirer une loi uniforme U U [0, 1], poser
Y = j si U rij ,
Xn +1 =
Xn = i sinon
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 10/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
Proposition (Réversibilité de Metropolis-Hastings) Si la
matrice de transition Q est irréductible et appériodique, alors la
chaîne (Xn ) produite par l’algorithme de Metropolis est
irréductible, appériodique et π est réversible pour cette chaîne.
Par ailleurs, tout ce qui a été dit dans le cas d’un espace d’états
…ni se généralise à un espace d’états in…ni, qu’il soit dénombrable
ou continu. Dans le cas d’un espace d’états continu, disons Rd ,
le but est de simuler selon la densité
fu ( x )
f (x ) = Cfu (x ) = R , (3)
fu (y ) dy
où, comme en formule 2, fu signi…e “f u unnormalized”.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 11/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
La matrice de transition Q = [qij ]1 i ,j M est remplacée par un
noyau de transition q (x, y ), c’est-à-dire que q (x, y ) 0 pour
tout couple (x, y ) et
Z
8 x 2 Rd , q (x, y ) dy = 1.
Comme son homologue discret, le noyau q peut se voir comme une
façon de se déplacer dans l’espace d’états Rd . La variable x étant
…xée, la quantité q (x, y ) s’interprète comme la densité de
probabilité d’aller en y sachant que l’on part du point x.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 12/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
Algorithme de Metropolis-Hastings, cas continu
Les hypothèses requises pour l’algorithme de Metropolis-Hastings
dans le cas continu sont alors :
- pour tout point x, on sait simuler selon la densité
q (x, ) ;
- pour tout couple (x, y ), on a q (x, y ) > 0 implique
q (y , x ) > 0 ;
- pour tout couple (x, y ) tel que q (x, y ) > 0, on sait
f (y )q (y ,x )
calculer f (x )q (x ,y ) .
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 13/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Généralisation de Metropolis-Hastings
Algorithme de Metropolis-Hastings, cas continu
Algorithme
1 partant de Xn = x, simuler Y = y selon la loi q (x, ) ;
2 calculer le rapport de Metropolis-Hastings
f (y )q (y , x )
r (x, y ) =
f (x )q (x, y )
3 tirer une loi uniforme U U [0, 1],
poser Xn +1 = Y = y si U r (x, y ) et Xn +1 = Xn = x sinon.
Sous des hypothèses raisonnables sur le noyau de transition q, la
chaîne (Xn ) a pour unique loi stationnaire f .
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 14/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Echantillonneur de Gibbs
Distribution de Boltzmann-Gibbs
Soit l’espace d’états S …ni (très grand), et supposons que l’on
veuille minimiser une fonction V : S ! R telle que V (x )
représente l’énergie associée à l’état x 2 S du système. Alors
sous des conditions assez générales, la distribution stationnaire
correspondant à l’équilibre thermodynamique des états du système
est une loi de Boltzmann-Gibbs :
1 V (x )
fT ( x ) = exp ,x 2 S (4)
ZT kT
dite aussi mesure de Gibbs associée au potentiel de V et à la
température T , où k représente la constante de Boltzmann et ZT
est la constante de normalisation, appelée fonction de répartition
en physique statistaique, elle est dé…nie par :
V (x )
ZT = ∑ exp kT
.
x 2S
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 15/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Echantillonneur de Gibbs
Distribution de Boltzmann-Gibbs
Nous supposons que la constante de Boltzmann a été incorporée
dans la température, ou ce qui revient au même, nous poserons
k = 1 dans la relation (4) ainsi que dans la fonction de répartition.
Alors l’équation devient
1 1
fT ( x ) = exp V (x ) , x 2 S. (5)
ZT T
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 16/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Echantillonneur de Gibbs
Lorsque la température tend vers zéro, la distribution alloue
donc toute la masse de probabilité aux états à énergie
minimale.
Puisque à basse température, la mesure de Gibbs se concentre
sur les minima globaux de V , une idée consiste à se donner
une suite de température (Tn ) décroissante vers 0 et pour
tout n, à simuler une chaîne de Markov (Xn ) qui a pour loi
stationnaire fT n .
Pour une matrice de transition Q (x, y ) donnée et véri…ant les
propriétés requises par Metropolis-Hastings (irréductibilité et
appériodicité), nous avons alors l’algorithme suivant, appelé
recuit cimulé (simulated annealing, venant de la
métallurgie).
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 17/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Echantillonneur de Gibbs
Algorithme de Metropolis pour simulated annealing
Algorithme
1 partant de Xn = x, simuler Y = y selon la loi Q (x, ) ;
2 calculer le rapport de Metropolis-Hastings
fT n ( y ) Q ( y , x )
rn (x, y ) =
fT n (x )Q (x, y )
Q (y , x ) 1
= exp (V (y ) V (x )) ,
Q (x, y ) Tn
3 tirer une loi uniforme U U [0, 1],
poser Xn +1 = Y = y si U rn (x, y ) et Xn +1 = Xn = x
sinon.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 18/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Echantillonneur de Gibbs
Remarque :
1 La chaîne n’est plus homogène en raison des changement de
températures Tn .
2 Si les transitions sont symétriques, i.e. Q (x, y ) = Q (y , x ) ,on
retrouve l’algorithme de Metropolis avec le rapport
1
rn (x, y ) = exp (E (y ) E (x )) ,
Tn
lequel est supérieur ou égal à 1 dès lors que E (y ) E (x ) .
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 19/20
Monté Carlo par Chaîne de Markov MCMC Algorithme de Metropolis-Hastings Echantillonneur de Gibbs
Echantillonneur de Gibbs
A…n de d’éviter de rester piégé dans un minimum local V , on
accepte systématiquement une transition vers une valeur de V plus
basse, mais on ne refuse pas systématiquement une transition vers
une valeur plus élevée.
On voit aussi que pour x et y …xés tels que V (y ) > V (x ) , plus la
température sera basse, moins on acceptera de transiter vers une
valeur élevée.
Chapitre 5 : Monté Carlo par Chaîne de Markov - Monte Carlo Markov Chain Dr. N. Arrar 20/20