Algorithme de Metropolis-Hastings

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}}$.

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é.

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.

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].

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).