Voir la source de Algorithme de Metropolis-Hastings

<!-- {{ébauche|statistiques}} -->
{{Infobox Méthode scientifique}} 

En [[statistique]], l''''algorithme de Metropolis-Hastings''' est une [[Méthode de Monte-Carlo par chaînes de Markov|méthode MCMC]] dont le but est d'obtenir un [[Échantillon (statistiques)|échantillonnage]] aléatoire d'une distribution de probabilité quand l'échantillonnage direct en est difficile. 

Plus formellement, étant donné une [[distribution de probabilité]] <math>\pi</math> sur un [[Univers (logique)|univers]] <math>\Omega</math>, cet [[algorithme]] définit une [[chaîne de Markov]] dont la [[Probabilité stationnaire d'une chaîne de Markov|distribution stationnaire]] est <math>\pi</math>. Il permet ainsi de tirer aléatoirement un élément de <math>\Omega</math> selon la loi <math>\pi</math>.

Un point essentiel de l'algorithme de Metropolis-Hastings est qu'il ne nécessite que la connaissance de <math>\pi</math> à une constante multiplicative près. En particulier, il n'est pas nécessaire de calculer la [[fonction de partition]] de <math>\pi</math>, 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éthode de rejet|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]]<ref>{{article
 | langue = en
 | nom1 = Metropolis | prénom1 = N.
 | nom2 =  Ulam | prénom2 = S.
 | année = 1949
 | titre = The Monte Carlo method
 | périodique = Journal of the American Statistical Association
 | volume = 44 | numéro = 247
 | pages = 335–341
 | doi = 10.2307/2280232 
}}.</ref> et sera décrite plus en détails en [[1953]] par Nicholas Metropolis, [[Arianna W. Rosenbluth]], [[Marshall Rosenbluth]], Augusta H. Teller, et [[Edward Teller]]<ref>{{article
 | langue = en
 | nom1 = Metropolis | prénom1 = N.
 | nom2 =  Rosenbluth | prénom2 = A.W.
 | nom3 = Rosenbluth | prénom3 = M.N.
 | nom4 =  Teller | prénom4 = A.H.
 | nom5 = Teller | prénom5 = E.
 | année = 1953
 | titre = Equations of State Calculations by Fast Computing Machines
 | périodique = Journal of Chemical Physics
 | volume = 21 | numéro = 6
 |pages = 1087–1092
}}.</ref>.
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, [[:en:W._K._Hastings|W. K. Hastings]] (1930-2016) a étendu l'algorithme au cas d'une distribution quelconque<ref>{{article
 | nom1 = Hastings | prénom1 = W.K.
 | année = 1970
 | titre = Monte Carlo Sampling Methods Using Markov Chains and Their Applications
 | périodique = Biometrika
 | volume = 57 | numéro = 1
 | pages = 97–109
}}.</ref>.

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''<ref name=Rosenbluth>{{article|titre=Genesis of the Monte Carlo Algorithm for Statistical Mechanics|auteur=M.N. Rosenbluth |journal=[[AIP Conference Proceedings]] | volume=690 | pages=22–30 | année=2003 | doi=10.1063/1.1632112 
}}</ref>
à l'occasion d'une conférence au [[Laboratoire national de Los Alamos|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<ref name=Gubernatis>{{article|titre=Marshall Rosenbluth and the Metropolis Algorithm |auteur=J.E. Gubernatis |journal=[[Physics of Plasmas]] | volume=12| pages=057303| année=2005| doi=10.1063/1.1887186 | bibcode=2005PhPl...12e7303G |numéro=5
|url=https://zenodo.org/record/1231899 }}</ref>.
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<ref name="Teller">Teller, Edward. ''Memoirs: A Twentieth-Century Journey in Science and Politics''. [[Perseus Publishing]], 2001, p. 328</ref>.
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 [[John von Neumann|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<ref name=Barth>Rosenbluth, Marshall. [https://www.aip.org/history-programs/niels-bohr-library/oral-histories/28636-1 "Oral History Transcript"]. American Institute of Physics</ref>, 
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<ref name=Dyson>{{article|titre=Marshall N. Rosenbluth|auteur=F. Dyson |journal=[[Proceedings of the American Philosophical Society]] | volume=250 | pages=404 | année=2006 
}}</ref>:
{{Citation étrangère|langue=en|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: {{citation|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 <math>x</math>, 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é <math>\pi</math>. Nous sommes dans le cas où il n'est pas simple de générer directement une suite de vecteurs suivant cette distribution <math>\pi</math>. Par ailleurs, il n'est pas nécessaire de connaître cette distribution, il suffit de connaître une fonction <math>f</math> telle que <math>f(x)</math> est proportionnelle à <math>\pi(x)</math>. 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 <math>\pi</math>.
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 <math>f</math> pour l'échantillon actuel à la valeur de <math>f</math> 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 <math>\pi</math>.

== 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 <math>f</math> une fonction proportionnelle à la distribution-cible, c'est-à-dire à la distribution de probabilité recherchée <math>\pi</math>.

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

=== 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 <math>\alpha</math> indique à quel point le nouvel état est probable étant donné l'état actuel, et selon la distribution <math>\pi</math>.
Si l'on cherche à se déplacer vers un état plus probable que l'état actuel (si <math>\alpha \ge 1 \ge u</math>), 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é <math>\pi</math> est élevée.
Par conséquent, la marche tend à visiter préférentiellement les régions de l'espace des états où la densité <math>\pi</math> 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 <math>g</math>.  <math>g(x,y)=g(x|y)</math> désigne la probabilité d'aller d'un état donné <math>y</math> vers un autre état donné <math>x</math>. Un processus markovien admet une unique [[Probabilité stationnaire d'une chaîne de Markov|distribution stationnaire]] <math>\pi(x)</math> si la chaîne de Markov est [[Graphe d'une chaîne de Markov et classification des états#Lexique_:_graphes-chaînes_de_Markov|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<ref>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.</ref>.

Soit <math>f</math> une fonction proportionnelle à la distribution-cible, c'est-à-dire à la distribution de probabilité recherchée <math>\pi</math>.

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

Il est important de noter que dans le cas général, le choix de la distribution <math>g</math> et du nombre d'itérations à réaliser<ref>Raftery, Adrian E., and Steven Lewis. "How Many Iterations in the Gibbs Sampler?" ''In Bayesian Statistics 4''. 1992.</ref> 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<ref name="Newman_Barkema">{{ouvrage|titre=Monte Carlo Methods in Statistical Physics |nom1=Newman |prénom1=M. E. J. |nom2=Barkema |prénom2=G. T. |année= 1999 |éditeur=Oxford University Press |lieu=USA |isbn=978-0198517979 }}</ref>.

=== 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 <math>\pi(x)</math> (qui peut être en grandes dimensions), on choisit une densité instrumentale conditionnelle <math>g</math> à 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<ref name=":0">{{article|titre = Adaptive Rejection Sampling for Gibbs Sampling |jstor = 2347565 |journal = Journal of the Royal Statistical Society. Series C (Applied Statistics) |date = 1992-01-01 |pages = 337–348 |volume = 41 |numéro = 2 |doi = 10.2307/2347565 |prénom1 = W. R. |nom1 = Gilks |prénom2 = P. |nom2 = Wild}}</ref> 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<ref>{{ouvrage|titre=Bayesian data analysis |date=2004 |éditeur=Chapman & Hall / CRC |champ libre=Gelman, Andrew |isbn=978-1584883883 |numéro d'édition= 2nd |lieu=Boca Raton, Fla. |oclc=51991499}}</ref>, 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 [[statistique multivariée|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èle bayésien hiérarchique|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<ref name=":0" />, 
l'algorithme adaptatif d'acceptation-rejet de Metropolis<ref>{{article|titre = Adaptive Rejection Metropolis Sampling within Gibbs Sampling |jstor = 2986138 |journal = Journal of the Royal Statistical Society. Series C (Applied Statistics) |date = 1995-01-01 |pages = 455–472 |volume = 44 |numéro = 4 |doi = 10.2307/2986138 |prénom1 = W. R. |nom1 = Gilks |prénom2 = N. G. |nom2 = Best |lien auteur2= Nicky Best |prénom3 = K. K. C. |nom3 = Tan}}</ref>, 
l'algorithme Metropolis–Hastings à 1 dimension, etc).

== Notes et références ==
{{Références}}
{{Traduction/Référence|lang1=en|art1=Metropolis–Hastings algorithm}}

== 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]]
* [[Échantillon (statistiques)|Échantillonnage]]
* [[Algorithme de Metropolis-Hastings à sauts réversibles]]

{{Portail|probabilités|physique|informatique théorique}}

[[Catégorie:Algorithme|Metropolis-Hastings]]
[[Catégorie:Processus stochastique]]
[[Catégorie:Physique statistique]]
[[Catégorie:Algorithme probabiliste|Metropolis Hastings]]
[[Catégorie:Méthode de Monte-Carlo par chaînes de Markov]]