Algorithme de Metropolis-Hastings

Algorithme de Metropolis-Hastings

Type	Algorithme
Inventeurs	Nicholas Metropolis, Marshall Rosenbluth, Edward Teller, W. K. Hastings (en), Enrico Fermi, Stanisław Ulam
Date d'invention	1953
Nommé en référence à	Nicholas Metropolis, W. K. Hastings (en)
Aspect de	Méthode de Monte-Carlo par chaînes de Markov

modifier - modifier le code - modifier Wikidata

En statistique, l'algorithme de Metropolis-Hastings est une méthode MCMC dont le but est d'obtenir un échantillonnage aléatoire d'une distribution de probabilité quand l'échantillonnage direct en est difficile.

Plus formellement, étant donné une distribution de probabilité ${\displaystyle \pi }$ sur un univers ${\displaystyle \Omega }$, cet algorithme définit une chaîne de Markov dont la distribution stationnaire est ${\displaystyle \pi }$. Il permet ainsi de tirer aléatoirement un élément de ${\displaystyle \Omega }$ selon la loi ${\displaystyle \pi }$.

Un point essentiel de l'algorithme de Metropolis-Hastings est qu'il ne nécessite que la connaissance de ${\displaystyle \pi }$ à une constante multiplicative près. En particulier, il n'est pas nécessaire de calculer la fonction de partition de ${\displaystyle \pi }$, tâche souvent difficile. Pour cette raison, cette méthode est très utilisée en physique statistique.

On peut noter que l'algorithme de Metropolis–Hastings (comme d'autres méthodes MCMC) est généralement utilisé pour l'échantillonnage de distributions multi-dimensionnelles, en particulier lorsque le nombre de dimensions est élevé. Pour les distributions unidimensionnelles, il existe habituellement d'autres méthodes pour générer des échantillons indépendants (par exemple les méthodes de rejet) qui permettent d'éviter les corrélations entre échantillons générés, problème inhérent aux méthodes MCMC.

Historique

Une première version de l'algorithme a été publiée en 1949 dans un article de Nicholas Metropolis et Stanisław Ulam[1] et sera décrite plus en détails en 1953 par Nicholas Metropolis, Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller, et Edward Teller[2]. Cette première version considérait le cas particulier de la distribution de Boltzmann, une des distributions les plus utilisées en physique statistique. En 1970, W. K. Hastings (1930-2016) a étendu l'algorithme au cas d'une distribution quelconque[3].

Il y a controverse sur la question de la paternité de l'algorithme. Metropolis connaissait les aspects computationnels de la méthode et avait introduit le terme Monte Carlo dans son article de 1949 avec Stanisław Ulam. Il a d'ailleurs dirigé le groupe qui a conçu l'ordinateur MANIAC I utilisé pour les expériences de 1952. Cependant, il n'existait pas de compte-rendu précis du développement de l'algorithme avant 2003. Peu avant sa mort, Marshall Rosenbluth décrivit le développement de l'algorithme au cours d'une présentation intitulée Genesis of the Monte Carlo Algorithm for Statistical Mechanics[4] à l'occasion d'une conférence au LANL qui célébrait le 50ème anniversaire de la publication de 1953. Gubernatis apporta de nouvelles clarifications dans un article de 2005 relatant cette conférence[5]. Rosenbluth y déclare qu'il avait réalisé le travail avec sa femme Arianna, et que Metropolis n'avait pas participé au développement de l'algorithme excepté pour leur accorder du temps de calcul.

Cela contredit le récit d'Edward Teller, qui relate dans ses mémoires que les cinq auteurs de l'article de 1953 y ont travaillé jour et nuit ensemble[6]. Le récit de Rosenbluth crédite Teller d'une suggestion cruciale: tirer avantage des moyennes d'ensemble utilisées en mécanique statistique plutôt que suivre une cinématique détaillée. C'est cette suggestion, dit Rosenbluth, qui l'a amené à réfléchir à une approche Monte Carlo généralisée, un sujet dont il aurait beaucoup discuté avec Von Neumann. Arianna Rosenbluth raconte à Gubernatis en 2003 qu'Augusta Teller avait commencé le code, mais qu'elle l'avait repris depuis le début. Dans un récit oral enregistré peu de temps avant sa mort[7], Rosenbluth déclare à nouveau que Teller avait proposé l'idée initiale, que lui (Rosenbluth) a résolu le problème et que sa femme Arianna l'a programmé.

Les compétences scientifiques reconnues de Rosenbluth apportent un certain crédit à son récit. Ainsi, Freeman Dyson écrit dans un texte biographique sur Rosenbluth[8]: « Many times I came to Rosenbluth, asking him a question [...] and receiving an answer in two minutes. Then it would usually take me a week of hard work to understand in detail why Rosenbluth's answer was right. He had an amazing ability to see through a complicated physical situation and reach the right answer by physical arguments. Enrico Fermi was the only other physicist I have known who was equal to Rosenbluth in his intuitive grasp of physics. » ce que l'on peut traduire par: « Souvent je posais une question à Rosenbluth à laquelle il répondait en deux minutes. Puis il me fallait une semaine de dur travail pour comprendre pourquoi la réponse de Rosenbluth était juste. Il avait une capacité extraordinaire à comprendre une situation physique complexe et à trouver la bonne réponse par des arguments physiques. Le seul autre physicien que j'ai connu qui avait une intuition physique équivalente était Enrico Fermi. »

Description intuitive

Nous voulons obtenir des tirages aléatoires d'un vecteur ${\displaystyle x}$, ayant un grand nombre de dimensions (pour une distribution à une dimension, on utilise d'autres algorithmes plus directs comme la méthode de rejet) selon une distribution de probabilité ${\displaystyle \pi }$. Nous sommes dans le cas où il n'est pas simple de générer directement une suite de vecteurs suivant cette distribution ${\displaystyle \pi }$. Par ailleurs, il n'est pas nécessaire de connaître cette distribution, il suffit de connaître une fonction ${\displaystyle f}$ telle que ${\displaystyle f(x)}$ est proportionnelle à ${\displaystyle \pi (x)}$. Cette propriété est particulièrement utile, car en pratique, il est souvent difficile de déterminer le facteur de normalisation.

L'algorithme de Metropolis–Hastings génère des valeurs d'échantillon de manière que plus on produit de valeurs, plus la distribution de ces valeurs se rapproche de la distribution recherchée ${\displaystyle \pi }$. Ces valeurs sont produites de manière itérative, la probabilité de l'échantillon à l'étape suivante ne dépend que de l'échantillon à l'étape considérée, la suite d'échantillons est donc une chaîne de Markov. Plus précisément, à chaque itération, l'algorithme choisit un candidat pour le prochain échantillon qui n'est basé que sur l'échantillon courant. Ensuite, le candidat est soit accepté avec une certaine probabilité (et dans ce cas il sera utilisé à l'étape suivante), soit rejeté avec la probabilité complémentaire (et dans ce cas l'échantillon actuel sera réutilisé à l'étape suivante). La probabilité d'acceptation est déterminée en comparant les valeurs de la fonction ${\displaystyle f}$ pour l'échantillon actuel à la valeur de ${\displaystyle f}$ pour l'échantillon candidat. Après un « grand nombre » d'itérations, la valeur des échantillons « perdent la mémoire » de l'état initial et suivent la distribution ${\displaystyle \pi }$.

Algorithme de Metropolis

On considère d'abord le cas de l'algorithme de Metropolis, qui est un cas particulier de l'algorithme de Metropolis–Hastings, où la probabilité de transition est symétrique. Hastings a généralisé cet algorithme à une probabilité de transition dissymétrique.

Pseudo-code

Soit ${\displaystyle f}$ une fonction proportionnelle à la distribution-cible, c'est-à-dire à la distribution de probabilité recherchée ${\displaystyle \pi }$.

1. Initialisation: choisir un point arbitraire ${\displaystyle x_{0}}$ pour être le premier échantillon et une probabilité de transition ${\displaystyle g}$ qui suggère un échantillon candidat ${\displaystyle x'}$ pour la valeur de l'échantillon suivant, connaissant la valeur de l'échantillon précédent ${\displaystyle x}$. On suppose ici que cette distribution de probabilité est symétrique, c'est-à-dire que ${\displaystyle g}$ doit vérifier ${\displaystyle g(x'\mid x)=g(x\mid x')}$. Un choix usuel pour une telle distribution est la distribution gaussienne centrée sur ${\displaystyle x}$, de manière que les points proches de ${\displaystyle x}$ aient une plus grande probabilité d'être visités à l'itération suivante[9].
2. A chaque itération t (échantillon courant ${\displaystyle x_{t}}$):
   • tirer aléatoirement un candidat ${\displaystyle x'}$ pour l'échantillon suivant selon la distribution ${\displaystyle g(x'\mid x_{t})}$.
   • Calculer le taux d'acceptation ${\displaystyle \alpha ={\frac {f(x')}{f(x_{t})}}}$, qui sera utilisé pour déterminer si le candidat est accepté ou rejeté. Comme ${\displaystyle f}$ est proportionnelle à la densité ${\displaystyle \pi }$, on a ${\displaystyle \alpha ={\frac {f(x')}{f(x_{t})}}={\frac {\pi (x')}{\pi (x_{t})}}}$.
   • Accepter ou rejeter:
     • Tirer un nombre uniformément aléatoire ${\displaystyle u\in [0,1]}$.
     • Si ${\displaystyle u\leq \alpha }$, alors accepter le candidat en effectuant ${\displaystyle x_{t+1}=x'}$ (donc si ${\displaystyle \alpha \geq 1}$, on accepte nécessairement),
     • Si ${\displaystyle u>\alpha }$, alors rejeter le candidat en effectuant ${\displaystyle x_{t+1}=x_{t}}$.

Interprétation

L'algorithme peut être interprété intuitivement de la manière suivante : à chaque itération, on tente de se déplacer dans l'espace des états possibles, le déplacement peut être accepté ou rejeté. Le taux d'acceptation ${\displaystyle \alpha }$ indique à quel point le nouvel état est probable étant donné l'état actuel, et selon la distribution ${\displaystyle \pi }$. Si l'on cherche à se déplacer vers un état plus probable que l'état actuel (si ${\displaystyle \alpha \geq 1\geq u}$), le déplacement est toujours accepté. Cependant, si l'on cherche à se déplacer vers un état moins probable que l'état actuel, alors le déplacement peut être rejeté et le rejet est d'autant plus probable que la chute de densité de probabilité ${\displaystyle \pi }$ est élevée. Par conséquent, la marche tend à visiter préférentiellement les régions de l'espace des états où la densité ${\displaystyle \pi }$ est élevée mais visite occasionnellement des régions de moindre densité.

Cas général

Généralités sur les processus de Markov

Un processus markovien est défini de manière unique par la distribution des probabilités de transition ${\displaystyle g}$. ${\displaystyle g(x,y)=g(x|y)}$ désigne la probabilité d'aller d'un état donné ${\displaystyle y}$ vers un autre état donné ${\displaystyle x}$. Un processus markovien admet une unique distribution stationnaire ${\displaystyle \pi (x)}$ si la chaîne de Markov est ergodique (ou irréductible, c'est-à-dire que tout état de l'espace est accessible depuis tout autre état), apériodique (le système ne revient pas dans le même état à intervalles réguliers) et récurrente positive (le nombre de pas pour revenir dans un même état est fini). On doit vérifier ces conditions pour pouvoir appliquer l'algorithme de Metropolis–Hastings.

Pseudo-code

Partiellement adapté de[10].

Soit ${\displaystyle f}$ une fonction proportionnelle à la distribution-cible, c'est-à-dire à la distribution de probabilité recherchée ${\displaystyle \pi }$.

1. Initialisation: choisir un point arbitraire ${\displaystyle x_{0}}$ pour être le premier échantillon et une probabilité de transition ${\displaystyle g}$ qui suggère un échantillon candidat ${\displaystyle x}$ pour la valeur de l'échantillon suivant, connaissant la valeur de l'échantillon précédent ${\displaystyle y}$. La condition de réversibilité de la chaîne de Markov implique que ${\displaystyle g(x\mid y)\pi (y)=g(y\mid x)\pi (x)}$.
2. A chaque itération t (échantillon courant ${\displaystyle x_{t}}$):
   • tirer aléatoirement un candidat ${\displaystyle x'}$ pour l'échantillon suivant selon la distribution ${\displaystyle g(x',x_{t})=g(x'\mid x_{t})}$.
   • Calculer le taux d'acceptation ${\displaystyle \alpha ={\frac {f(x')g(x_{t},x')}{f(x_{t})g(x',x_{t})}}}$, qui sera utilisé pour déterminer si le candidat est accepté ou rejeté. Comme ${\displaystyle f}$ est proportionnelle à la densité ${\displaystyle \pi }$, on a ${\displaystyle \alpha ={\frac {f(x')g(x_{t},x')}{f(x_{t})g(x',x_{t})}}={\frac {\pi (x')g(x_{t},x')}{\pi (x_{t})g(x',x_{t})}}}$.
   • Accepter ou rejeter:
     • Tirer un nombre uniformément aléatoire ${\displaystyle u\in [0,1]}$.
     • Si ${\displaystyle u\leq \alpha }$, alors accepter le candidat en effectuant ${\displaystyle x_{t+1}=x'}$ (donc si ${\displaystyle \alpha \geq 1}$, on accepte nécessairement),
     • Si ${\displaystyle u>\alpha }$, alors rejeter le candidat en effectuant ${\displaystyle x_{t+1}=x_{t}}$.

Il est important de noter que dans le cas général, le choix de la distribution ${\displaystyle g}$ et du nombre d'itérations à réaliser[11] sont des paramètres libres de la méthode et doivent être adaptés au problème considéré. Dans le cas d'un espace discret, le nombre d'itérations doit être de l'ordre du temps d'autocorrélation du processus de Markov[12].

Analyse et comparaison à d'autres méthodes

L'algorithme de Metropolis–Hastings est une des méthodes les plus générales dans la famille des méthodes MCMC, dans le sens où il impose très peu de conditions sur la densité cible. À partir de la densité cible ${\displaystyle \pi (x)}$ (qui peut être en grandes dimensions), on choisit une densité instrumentale conditionnelle ${\displaystyle g}$ à partir de laquelle il est relativement facile de simuler.

Par comparaison avec d'autres algorithmes d'échantillonnage tel que la méthode de rejet et sa variante adaptative[13] qui génère directement des échantillons indépendants depuis la distribution, l'algorithme de Metropolis–Hastings (et les méthodes MCMC en général) présentent des désavantages:

• Les échantillons sont corrélés. Un ensemble de points proches au cours du processus seront corrélés entre eux et ne reflète donc pas correctement la distribution. Par conséquent, si l'on souhaite avoir des échantillons indépendants, on doit éliminer la majorité des échantillons pour ne conserver qu'un échantillon tous les n pas, avec n « suffisamment grand ». On peut déterminer n par exemple en mesurant l'autocorrélation entre échantillons proches. L'autocorrélation peut être réduite en augmentant la largeur des sauts entre deux échantillons, mais augmenter la taille des sauts augmente aussi la probabilité de rejet du nouvel échantillon. Dans le cas où le processus nécessite un grand nombre de sauts pour éliminer l'autocorrélation, on dit que la chaîne de Markov a un long temps de mélange.
• Bien que la chaîne de Markov converge vers la distribution souhaitée, les échantillons initiaux peuvent suivre une distribution très différente, en particulier s'ils se trouvent dans une région de faible densité de probabilité. Par conséquent, une période de rodage est souvent nécessaire[14], et les premiers échantillons seront éliminés.

En revanche, les méthodes de rejet les plus simples peuvent souffrir du « fléau de la dimension », qui se traduit ici par le fait que la probabilité de rejet augmente exponentiellement avec le nombre de dimensions. Les méthodes MCMC n'ont pas ce problème jusqu'à un certain point et sont donc souvent les seules solutions disponibles lorsque le nombre de dimensions de la distribution à échantillonner est élevé.

Dans le cas des distributions multivariées, l'algorithme de Metropolis–Hastings tel que décrit plus haut suppose de choisir à chaque itération un nouvel échantillon dans l'espace multi-dimensionnel. Quand le nombre de dimensions est élevé, trouver le pas du déplacement peut être difficile, car chaque dimension peut se comporter très différemment, or le pas du saut doit être adapté à chaque dimension pour éviter un mélange trop lent. Une approche alternative qui peut mieux fonctionner dans de telles situations est l'échantillonnage de Gibbs. Cette approche consiste à choisir un nouvel échantillon pour chaque dimension séparément des autres. Elle est utilisable en particulier lorsque la distribution multivariée est composée d'un ensemble de variables aléatoires conditionnées par un petit ensemble d'autres variables, ce qui est typiquement le cas dans les modèles bayésiens hiérarchiques. Les variables individuelles sont échantillonnées une à une et chaque variable est conditionnée sur la valeur la plus récente de celles dont elle dépend. Plusieurs algorithmes peuvent être utilisés pour choisir ces échantillons individuels, en fonction de la forme précise de la distribution multivariée (par exemple l'algorithme adaptatif d'acceptation-rejet[13], l'algorithme adaptatif d'acceptation-rejet de Metropolis[15], l'algorithme Metropolis–Hastings à 1 dimension, etc).

Notes et références

1. (en) N. Metropolis et S. Ulam, « The Monte Carlo method », Journal of the American Statistical Association, vol. 44, n^o 247,‎ 1949, p. 335–341 (DOI 10.2307/2280232).
2. (en) N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller et E. Teller, « Equations of State Calculations by Fast Computing Machines », Journal of Chemical Physics, vol. 21, n^o 6,‎ 1953, p. 1087–1092.
3. W.K. Hastings, « Monte Carlo Sampling Methods Using Markov Chains and Their Applications », Biometrika, vol. 57, n^o 1,‎ 1970, p. 97–109.
4. M.N. Rosenbluth, « Genesis of the Monte Carlo Algorithm for Statistical Mechanics », AIP Conference Proceedings, vol. 690,‎ 2003, p. 22–30 (DOI 10.1063/1.1632112)
5. J.E. Gubernatis, « Marshall Rosenbluth and the Metropolis Algorithm », Physics of Plasmas, vol. 12, n^o 5,‎ 2005, p. 057303 (DOI 10.1063/1.1887186, Bibcode 2005PhPl...12e7303G, lire en ligne)
6. Teller, Edward. Memoirs: A Twentieth-Century Journey in Science and Politics. Perseus Publishing, 2001, p. 328
7. Rosenbluth, Marshall. "Oral History Transcript". American Institute of Physics
8. F. Dyson, « Marshall N. Rosenbluth », Proceedings of the American Philosophical Society, vol. 250,‎ 2006, p. 404
9. Dans l'article original, ${\displaystyle g(x'\mid x)}$ est la distribution de Boltzmann, car l'algorithme était utilisé dans le contexte de la mécanique statistique.
10. Vanessa Bergeron Laperrière (Été 2010), (supervisée par Mylène Bédard), L’Algorithme Metropolis–Hastings Projet de recherche CRSNG, Département de Mathématiques et Statistique Université de Montréal.
11. Raftery, Adrian E., and Steven Lewis. "How Many Iterations in the Gibbs Sampler?" In Bayesian Statistics 4. 1992.
12. M. E. J. Newman et G. T. Barkema, Monte Carlo Methods in Statistical Physics, USA, Oxford University Press, 1999 (ISBN 978-0198517979)
13. W. R. Gilks et P. Wild, « Adaptive Rejection Sampling for Gibbs Sampling », Journal of the Royal Statistical Society. Series C (Applied Statistics), vol. 41, n^o 2,‎ 1^er janvier 1992, p. 337–348 (DOI 10.2307/2347565, JSTOR 2347565)
14. Bayesian data analysis (Gelman, Andrew), Boca Raton, Fla., Chapman & Hall / CRC, 2004, 2nd éd. (ISBN 978-1584883883, OCLC 51991499)
15. W. R. Gilks, N. G. Best et K. K. C. Tan, « Adaptive Rejection Metropolis Sampling within Gibbs Sampling », Journal of the Royal Statistical Society. Series C (Applied Statistics), vol. 44, n^o 4,‎ 1^er janvier 1995, p. 455–472 (DOI 10.2307/2986138, JSTOR 2986138)

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Metropolis–Hastings algorithm » (voir la liste des auteurs).

Voir aussi

Articles connexes

• Chaîne de Markov
• Recuit simulé
• Optimisation combinatoire
• Méthode de Monte-Carlo par chaînes de Markov
• Méthode de Monte-Carlo
• Échantillonnage
• Algorithme de Metropolis-Hastings à sauts réversibles