Algorithmes MCMC

Algorithmes MCMC : principe général

Le principe général des algorithmes MCMC est le suivant : pour produire une approximation acceptable d’une intégrale ou d’une autre fonctionnelle d’une distribution d’intérêt (i.e. la loi a posteriori), il suffit de générer une chaîne de Markov dont la distribution limite est la distribution d’intérêt (i.e. la loi a posteriori), puis d’y appliquer la méthode de Monte-Carlo.

Il faut donc avoir une double convergence :

1. la convergence de la chaîne de Markov vers sa distribution stationnaire : \(\forall X_0,\, X_n\underset{{n\to+\infty}}{\xrightarrow{\quad\mathcal{L}\quad}} \tilde{p}\)
2. la convergence de la Loi des Grands Nombres, une fois la distribution stationnaire atteinte : \(\frac{1}{N}\sum_{i=1}^N f(X_{n+i}) \underset{N\to+\infty}{\xrightarrow{\qquad~~}} \mathbb{E}[f(X)]\)

\[\overbrace{X_0\rightarrow X_1 \rightarrow X_2\rightarrow \dots \rightarrow X_n}^\text{convergence de la chaîne de Markov}% \rightarrow \overbrace{X_{n+1} \rightarrow X_{n+2} \rightarrow \dots \rightarrow X_{n+N}}^\text{échantillon de Monte-Carlo}\] Cette double convergence est garantie par le théorème ergodique.

Les algorithmes MCMC utilisent une approche d’acceptation-rejet :

\(\qquad\)Schéma général des algorithmes MCMC

Pour la loi instrumentale de proposition \(q\) il n’existe pas de choix absolument optimal mais une infinité de lois possibles (certaines meilleures que d’autres). Afin de garantir la convergence vers la loi cible \(\tilde{p}\) : (i) le support de \(q\) doit contenir le support \(\tilde{p}\), (ii) \(q\) ne doit pas générer de valeurs périodiques. Idéalement, on choisit \(q\) de manière à ce que son calcul soit simple (par exemple on peut choisir \(q\) symétrique).

L’algorithme de Metropolis-Hastings

L’algorithme de Metropolis-Hastings est un algorithme très simple et très général permettant d’échantillonner selon des lois uni- ou multi-dimensionnelles.

On peut reformuler la probabilité d’acceptation \(\alpha^{(t)}\) ainsi : \(\displaystyle\alpha^{(t)} = \min\left\{1, \frac{\tilde{p}(y^{(t)})}{\tilde{p}(x^{(t-1)})}\frac{q(x^{(t-1)}|y^{(t)})^{(t)}}{q(y^{(t)}|x^{(t-1)})}\right\}\). On voit donc qu’on peut la calculer en ne connaissant \(\tilde{p}\) qu’à une constante près, puisqu’elle se simplifie dans ce ratio.

Dans certains cas particuliers (très utilisés en pratique), le calcul de \(\alpha^{(t)}\) est simplifié :

• Metropolis-Hastings indépendant : \(q(y^{(t)}|x^{(t-1)}) = q(y^{(t)})\)
• Metropolis-Hastings à marche aléatoire : \(q(y^{(t)}|x^{(t-1)}) = g(y^{(t)}-x^{(t-1)})\). Si \(g\) est symétrique (\(g(-x) = g(x)\)), alors le calcul de la probabilité d’acceptation \(\alpha^{(t)}\) se simplifie : \(\frac{\tilde{p}(y^{(t)})}{\tilde{p}(x^{(t-1)})}\frac{q(y^{(t)}|x^{(t-1)})}{q(x^{(t-1)}|y^{(t)})} = \frac{\tilde{p}(y^{(t)})}{\tilde{p}(x^{(t-1)})}\frac{\cancel{g(y^{(t)}-x^{(t-1)})}}{\cancel{g(x^{(t-1)}-y^{(t)})}} = \frac{\tilde{p}(y^{(t)})}{\tilde{p}(x^{(t-1)})}\)

L’algorithme de Metropolis-Hastings est un algorithme très simple et très général permettant d’échantillonner de manière uni- ou multi-dimensionnelle. Le choix de la distribution instrumentale est crucial, mais difficile, et a un impact considérable sur les performances de l’algorithme (un fort taux de rejet implique souvent des temps de calculs très importants). De plus, c’est un algorithme qui devient inefficace dans les problèmes de trop grande dimension. L’algorithme du recuit-simulé ainsi que l’échantillonneur de Gibbs sont des algorithmes permettant en partie de pallier à certaines de ces limites.

L’algorithme du recuit-simulé

Afin de dépasser certaines limitations de l’algorithme de Metropolis-Hastings, ici l’idée est de faire varier le calcul de la probabilité d’acceptation \(\alpha^{(t)}\) au cours de l’algorithme. La probabilité d’acceptation doit d’abord être grande afin de bien explorer l’ensemble de l’espace d’état, puis diminuer au fur et à mesure que l’algorithme converge vers une région de l’espace, afin que les nouvelles valeurs acceptées se concentrent autour du mode de convergence. Cela consiste à introduire dans l’algorithme de Métropolis-Hastings une “température” variant à chaque itération et notée \(T(t)\) :

Par exemple, on peut prendre \(T(t) = T_0 \left(\frac{T_f}{T_0}\right)^\frac{t}{n}\) avec \(T_0\) la température de base, \(n\) le nombre d’itérations au-delà duquel on pense être proche de la convergence, \(T_f\) la température après \(n\) itérations. Cet algorithme est particulièrement utile lors de la présence d’optimums locaux.

Échantillonneur de Gibbs

Lorsque la dimension (de \(x\)) augmente, il devient très difficile de proposer des valeurs probables dans les algorithmes utilisant la stratégie d’acceptation-rejet. L’idée de l’échantillonneur de Gibbs est de générer \(x\) coordonnée par coordonnée, en conditionnant sur les dernières valeurs obtenues. Il faut donc que \(x\) admette une décomposition telle que \(x=(x_1, \dots, x_d)\), et que les distributions \(p(x_i | x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d)\) soit connues et aisément simulables. L’échantillonneur de Gibbs est une suite d’étapes de Metropolis-Hastings (coordonnée par coordonnée) où les propositions échantillonnées sont toujours acceptée (\(\alpha = 1\)). On obtient cette acceptation inconditionnelle en imposant les lois de propositions : il s’agit de la distribution conditionnelle respective de chaque coordonnée. On peut donc voir l’échantillonneur de Gibbs comme un algorithme de réactualisation composante par composante :

Remarque : si l’on ne connaît pas certaines lois conditionnelles pour certaines coordonnées, on peut tout de les échantillonner en introduisant une étape d’acceptation-rejet pour cette coordonnée uniquement. On parle alors d’algorithme de Métropolis à l’intérieur de Gibbs (Metropolis within Gibbs).