\section{La méthode MCMC}

Ce cours-td est issu de quatres sources (\cite{BMP}, \cite{Br},\cite{Hag} et \cite{Lig}). Le meilleur exposé est sans nul doute \cite{Hag}. Néanmoins, pour ce qui est des motivations de l'étude du modèle d'Ising, \cite{Lig} est nécessaire même si \cite{Br} l'aborde d'une manière moins agréable. Ce qui concerne les mesures de Gibbs en volume infini est donné sans preuves (cf. \cite{Lig} pour de plus amples détails). La seule raison d'exposer ces résultats est de donner une motivation non dénuée d'intérêt à la simulation d'une mesure de Gibbs. D'autre part, la ``convergence'' de l'algorithme de Propp-Wilson n'est pas donné dans \cite{Hag} mais il l'est dans \cite{BMP}, exercice 4.38 p.123.


\subsection{Les méthodes classiques}

On peut distinguer trois grandes méthodes de simulations de variables aléatoires.\\
1. Méthode de l'inverse : On cherche à simuler une variable aléatoire $Z$ de fonction de répartition $F$. Pour ce faire, on simule une variable aléatoire uniforme sur $[0,1]$. Soit $F^{-1}$ la pseudo inverse de $F$ alors $F^{-1}(U)$ a pour loi $F$.\\
2. Méthode de rejet : On cherche à simuler une variable aléatoire $Z$ ayant une densité $f$. On suppose que l'on est capable de générere une suite i.i.d. de variables aléatoires $(Y_n)_n$ ayant une loi à densité $g$ telle que
$$\cfrac{f(x)}{g(x)}\leq c$$
Soit $(U_n)_n$ une suite i.i.d. de variables aléatoires uniformément distribuées sur $[0,1]$. On définit $\tau$ comme le premier instant $n$ pour lequel on a $$U_n \leq \cfrac{f(Y_n)}{cg(Y_n)}$$
et
$$Z=Y_{\tau}$$
Alors $Z$ est une loi ayant pour densité $f$.\\
3. Méthode MCMC (Monte Carlo Markov Chain) : c'est celle que nous allond décrire plus en détail maintenant.\\

Avant de décrire la méthode MCMC (dans le cas d'une loi discrète), commençons par donner quelques critiques en ce qui concerne les deux premières méthodes. Les méthodes de l'inverse et de rejet ont des contre-parties quand $Z$ est une variable discrète de loi $\pi$ sur un ensemble fini $E$ supposé très grand.\\
La méthode de l'inverse consiste à générer une variable uniforme $U$ sur $[0,1]$ et à poser $Z=i$ ssi $\sum_{l=1}^{i-1} \pi_l \leq U <\sum_{l=1}^{i} \pi_l$. Quand la taille de $E$ est très grande, on a une partition de l'intervalle $[0,1]$ constituée de petits intervalles. Savoir où $U$ se trouve exactement demande donc un grand coût en précision de calculs. Dans une simulation où $E$ n'est pas une suite d'entiers mais un espace de configurations plus ou moins compliqué, il est nécessaire de définir un codage de l'espace des configurations (et de coder et décoder) ce qui peut être assez lourd.\\
La méthode de rejet dans le cas discret utilise une distribution $p$ sur $E$ telle que $\cfrac{\pi_i}{p_i}\leq c$ pour tout $i\in E$. Le problème vient du fait que bien souvent $\pi$ est connue à une constante de normalisation près si bien que la méthode de rejet (tout comme la méthode de l'inverse) devient inutilisable.\\

\subsection{Le principe de la méthode MCMC}

La méthode MCMC est la suivante. On construit une chaîne de Markov irréductible aperiodique sur $E$ avec $\pi$ pour distribution stationnaire. Puisque $E$ est fini, la chaîne est ergodique et pour tout $i\in E$ et toute distribution initiale $\mu$,
$$\lim_{n\rightarrow \infty} \PP_{\mu} (X_n =i)=\pi_i$$
Aussi, quand $n$ est ``grand, on peut considérer que $X_n$ a une distribution proche de $\pi$. Bien entendu, il se pose deux problème majeur dans cette méthode :\\
1) Comment construire une chaîne $X$ ayant $\pi$ pour distribution stationnaire ?\\
2) Comment évaluer si la distribution de $X_n$ est proche de $pi$ ? Ou plus exactement, comment obtenir une estimation de la vitesse de convergence de la loi de $X_n$ vers $\pi$ en fonction de $n$ ?\\

Un élément de réponse à la première question est fournie par les algorithme de Métropolis (1953) et Barker (1965). On remarquera au passage que ce dernier algorithme ne nécessite la connaisance de $\pi$ qu'à un facteur multiplicatif près. Pour ce qui est de la deuxième question, on sait que la vitesse de convergence est géométrique et donnée par la deuxième plus grande valeur propre de la matrice de transition. Tout revient alors à obtenir une bonne estimation de cette seconde valeur propre... On trouvera tous les compléments utiles dans le chapitre $6$ de $\cite{Br}$.

\dexer
Algorithme de Metropolis ($\cite{BMP}$ exercice 4.37)\\
Soit $E$ un ensemble d'\'etats fini et $P$ une matrice de transition sur $E$ sym\'etrique et irr\'eductible. Soit $\pi$ une probabilit\'e sur $E$ telle que $\pi(x)>0$ pour tout $x\in E$. On d\'efinit sur $E$ une nouvelle matrice de transition $Q$ par
\begin{equation*}
Q(x,y)=
\begin{cases}
P(x,y) \text {   si $x\neq y$ et $\pi(y)\geq \pi (x)$}\\
P(x,y)\cfrac{\pi(x)}{\pi(y)} \text{   si $\pi(y)< \pi(x)$}\\
1-\sum_{z\neq x} Q(x,z) \text{   si $y=x$}
\end{cases}
\end{equation*}
\quest a) Montrer que la loi $\pi$ est réversible pour $Q$ et donc invariante.

\quest b) Montrer que $Q$ est irréductible.

\quest c) On suppose que $\pi$ n'est pas la distribution uniforme (sinon on aurait $Q=P$). On note $M$ l'ensemble des états $x \in E$ tels que $\pi (x)= \sup_{y \in E} \pi (y)$.

\quest c1) Montrer qu'il existe $x_0 \in M$ tel que $P(x_0, y)>0$ pour un certain $y\notin M$. En déduire que $Q (x_0, z) < P(x_0, z)$ pour un $z \in E$ au moins. Monter que $Q(x_0, x_0)>0$.

\quest c2) Montrer que $Q$ est apériodique.


On a donc montré que $\pi$ est r\'eversible pour $Q$ donc invariante et que $Q$ est irr\'eductible et aperiodique. Par le th\'eor\`eme ergodique, on est assur\'e que pour $n$ assez grand, la loi de la chaîne de Markov $(X_n)_n$ de matrice de transition $Q$ est proche de $\pi$. Cette m\'ethode de simulation est particuli\`erement int\'eressante lorsque l'on ne conna\^ {i}t $\pi$ qu'\`a une constante multiplicative pr\`es (c'est ce qui arrive fr\'equemment en m\'ecanique statistique o\`u la fonction de r\'epartition est inconnue, ou du moins si $E$ est si grand que la calculer a un co\^ut trop grand). Une exellente application, malheureusement en dehors des thèmes du programme et assez difficile mais faisable, de l'algoritme de Metropolis est l'illustration de la transition de phase dans le modèle d'Ising (Cf. $\cite{Br}, \cite{Lig}$).\\
\fexer


\subsection{Simulation d'une chaîne de Markov}

La question est la suivante. Comment simuler une chaîne de Markov $(X_0,\ldots,X_n,\ldots)$ à espace d'états fini $S=\{s_1,\ldots,s_k\}$, de matrice de transition $P$, avec la distribution $\mu$ pour loi initiale. Nous admettons que nous pouvons simuler une suite $(U_0,\ldots,U_n,\ldots)$ (en fait, une suite finie...) de variables aléatoires uniformes indépendantes sur $[0,1]$. Les deux ingrédients principaux sont la fonction d'initiation $\psi$ et la fonction de mise à jour $\phi$. On les définit de la manière suivante.
\begin{equation}
\psi(x)=
\begin{cases}
s_1 \qquad \mbox{pour $x\in [0, \mu (s_1))$}\\
s_2 \qquad \mbox{pour $x\in [\mu (s_1), \mu (s_1)+\mu (s_2))$}\\
\cdot\\
\cdot\\
s_i \qquad \mbox{pour $x\in [\sum_{j=1}^{i-1} \mu (s_j), \sum_{j=1}^i \mu (s_j))$}\\
s_k \qquad \mbox{pour $x\in [\sum_{j=1}^{k-1} \mu (s_j),1]$}
\end{cases}
\end{equation}
et pour $i\in \{1,\ldots,k\}$,
\begin{equation}
\phi(s_i,x)=
\begin{cases}
s_1 \qquad \mbox{pour $x\in [0, P_{i,1})$}\\
s_2 \qquad \mbox{pour $x\in [P_{i,1}, P_{i,1}+P_{i,2})$}\\
\cdot\\
\cdot\\
s_j \qquad \mbox{pour $x\in [\sum_{l=1}^{j-1} P_{i,l}, \sum_{l=1}^j P_{i,l})$}\\
s_k \qquad \mbox{pour $x\in [\sum_{l=1}^{k-1} P_{i,l},1]$}
\end{cases}
\end{equation}

Soit $(U_0,\ldots,U_n,\ldots)$ une suite de variables aléatoires uniformes sur $ [0,1]$ indépendantes. On pose $X_0=\psi (U_0), X_1 = \phi (X_0, U_1),\ldots, X_n=\phi (X_{n-1}, U_n), \ldots)$.\\

\quest 1) Montrer que $(X_n)_n$ est une chaîne de Markov de matrice de transition $P$, de loi initiale $\mu$.
\quest 2) En déduire un algorithme pour simuler une chaîne de Markov donnée.

\subsection{Un modèle d'exclusion}

Soit $\Lambda_{m} =\{1,\ldots,m\}^2$ un carré dans $\ZZ^2$ de longuer $m$. Une configuration $\zeta$ sur $\Lambda_m$ est une fonction de $\Lambda_m$ dans $\{0,1\}$. Une configuration $\zeta$ est dite faisable si pour tout $x$ et $y$ de $\Lambda_m$, $|x-y|=1$ implique $\zeta(x)\eta (y) =0$ (autrement dit, on ne peut pas avoir deux sites adjacents occupés avec des $1$). Soit $Z_m$ le nombre de configurations faisables. On définit une probabilité $\mu_m$ sur l'espace $\Omega_m=\{0,1\}^m$ des configurations par:
\begin{equation}
\mu_m (\zeta)=Z_m^{-1} \indic_{\{\zeta \; \mbox{faisable}\}}
\end{equation}


\quest 1) Soit $ n(\zeta)$ le nombre de $1$ dans la configuration $\zeta$. On s'intéresse à l'estimation de $\mu (n)$. Pour cela, on doit simuler une variable aléatoire $X$ ayant pour loi (ou approximativement) la loi $\mu$. Admettons que l'on dispose d'une procédure permettant de simuler une telle variables aléatoire. Comment faire pour estimer (en terme de simulation) la quantité $\mu (n)$ ?\\
\quest 2) On se propose maintenant de donner un algorithme MCMC pour ce modèle d'exclusion, i.e. simuler par la méthode MCMC une variable aléatoire $X$ ayant pour loi $\mu_m$.\\
\subquest a) Pourquoi une méthode MCMC s'avère nécessaire pour $m$ grand ?\\
\subquest b) Vérifier que l'algorithme suivant convient:\\
1. Choisir un site $s$ dans $\Lambda_m$ de manière aléatoire et uniforme.\\
2. Simuler une variable aléatoire de Bernoulli $U$ de paramètre $1/2$.\\
3. Si $U=1$ et si tous les voisins de $s$ ont pour valeur $0$ alors mettre la nouvelle valeur de $s$ à $1$, sinon mettre à $0$.\\
4. Pour les autres sites (i.e. différents de $s$), ne rien faire.\\
\subquest c) Le simuler.\\
\quest 3) Donner une interprétation en termes de ``modèle de télécommunication'' et de ``dynamique des populations'' de cette mesure.

\section{Critique de la méthode MCMC}

On peut adresser deux grandes critiques à la méthodes MCMC.\\
A) L'argument théorique qui fonde la méthode MCMC est le théorème qui dit que la loi $\mu_n$ au temps $n$ d'une chaîne de Markov irréductible aperiodique (partant de n'importe quel état) converge vers sa distribution stationnaire $\pi$ (distribution que l'on cherche à simuler). Bien entendu, ce n'est qu'une convergence et cela n'implique pas l'existence d'un $n$ pour lequel $\mu_n =\pi$. Il y a donc toujours une erreur lorsque l'on approxime $\pi$ par $\mu_n$.\\
B) Pour rendre cette erreur petite, on doit être capable d'une manière ou d'une autre d'estimer cette erreur (erreur que l'on mesure au moyen de la distance en variation).\\
Le problème soulevé en A) est un problème inhérent à toute simulation numérique et on s'en accomode aisément, d'autant plus que l'ordinateur lui-même ne donne que des résultats approximatifs (quelque soit la précision dont il dispose, $\RR$ lui reste terra incognita). Le second problème B) est quant à lui beaucoup plus serieux. Il s'est avérer qu'estimer cette erreur (i.e. majorer $d(\mu_n, \pi)$) était très difficile. Nous disposons cependant de résultats dans ce sens (théorème $8.1$ p.56
de $\cite{Hag}$ par exemple). En pratique, on opère en général de la manière suivante. On simule la chaîne de Markov jusqu'à un temps $n$ de l'ordre $10^5$ et on croise les doigts! Dans la section suivante, nous présentons un algorithme ``parfait'' dans le sens où il simule une loi qui est exactement égale à la distribution $\pi$ cherchée. Bien entendu, cet algorithme possède lui aussi des désavantages que nous expliciterons.


\section{L'algorithme de Propp-Wilson}

On adpote les notations classsiques de ce chapitre. Soit une chaîne de Markov $(X_0,\ldots,X_n,\ldots)$ à espace d'états fini $S=\{s_1,\ldots,s_k\}$, de matrice de transition $P$, avec la distribution $\mu$ pour loi initiale. La fonction de mise à jour est notée $\phi:S\times [0,1] \rightarrow S$.\\
On considère une suite de nombre $N_1=1,\ldots, N_m=2^{m-1},\ldots$. Soit $(U_0, U_{-1}, \ldots,U_{-n},\ldots)$ une suite de variables aléatoires uniformes indépendantes sur $[0,1]$. L'algorithme est le suivant:\\
\noindent
1. Prendre $m=1$.\\
2. Pour tous les $s \in \{s_1, \ldots, s_k\}$, simuler la chaîne de Markov partant de $-N_m$ dans l'état $s$ jusqu'au temps $0$ en utilisant $\phi$ et les variables aléatoires $U_{-N_m +1}, \ldots, U_{-1}, U_{0}$ (ce sont les mêmes pour les $k$ trajectoires simulées).\\
3. Si les $k$ chaînes dans l'étape $2$ se terminent dans le même état $s'$ au temps $0$ alors l'algorithme se termine et donne $s'$. Autrement, il passe à l'étape $4$.\\
4. Incrémenter $m$ de $1$ et continuer avec l'étape $2$.\\

Quelques commentaires sur cet algorithme. Première chose, il est important de dire que à la $m$-ème étape, l'algorithme réutilise certaine les variables aléatoires $(U_{-N_{m-1} +1}, \ldots, U_0)$ de l'étape $m-1$. Ceci est économique en appelle de simulation de variables aléatoires uniformes mais coûteux en mémoire puisque l'on doit conserver en mémoire les variables aléatoires du passé. Néanmoins, ceci est nécessaire pour que l'algorithme marche (un exemple est donné dans \cite{Hag} où l'algorithme ne fonctionne pas si l'on ne réutilise pas les anciennes uniformes). Deuxièmement, a priori, l'algorithme contient une boucle infinie. Il peut très bien ne jamais terminer. En fait, un résultat (que nous ne prouverons pas mais que l'on peut trouver dans \cite{BMP}, exercice 4.38) dit que l'on peut choisir une bonne fonction $\phi$, que nous dirons \textit{valide},  telle que avec probabilité $1$, l'algorithme termine toujours. Mais un mauvais choix de la fonction $\phi$ peut conduire à une boucle infinie (cf. le problème élémentaire 10.2 de $\cite{Hag}$). Troisièmement, il peut sembler étrange de partir dans le ``passé'', en plus à des temps en $-2^m$. Boh! Nous passons maintenant à la preuve du théorème suivant:

\begin{thf}
Soit $P$ la matrice de transition d'une chaîne de Markov irréductible et aperiodique à espace d'états $S=\{s_1,\ldots, s_k\}$ et ayant pour distribution stationnaire $\pi=(\pi_1,\ldots,\pi_k)$. Soit $\phi$ une fonction de mise à jour \textit{valide} pour $P$ et considérons l'algorithme de Propp-Wilson. On suppose que l'algorithme se termine avec probabilité $1$ et on note $Y$ la sortie de l'algorithme. Alors pour tout $i\in \{1,\ldots, k\}$,
$$\PP (Y=s_i) = \pi_i$$
\end{thf}

\quest 1) Fixons $s_i \in S$. Il suffit de prouver que pour tout $\epsilon >0$, on a
$$| \PP(Y =s_i) - \pi_i| \leq \epsilon$$
Montrer qu'il existe $M$ tel que
$$\PP(\mbox{l'algorithme n'a pas besoin de faire plus de M étapes}) \geq 1-\epsilon$$

\quest 2) On fixe un tel $M$. On considère une chaîne de Markov $\tilde{Y}$ partant du temps $-N_M$ avec la même fonction de mise à jour $\phi$ et la même suite d'uniformes $U_{-N_M +1}, \ldots, U_0$ mais avec pour loi initiale (au temps $-N_{M}$) la probabilité stationnaire $\pi$. Monter que
$$\PP (Y \neq {\tilde{Y}}_0) \leq \epsilon$$

\quest 3) En déduire $| \PP(Y =s_i) - \pi_i| \leq \epsilon$.\\


Passons maintenant à la critique de cet algorithme. Nous nous intéressons à la simulation de distributions $\pi$ sur $S$ lorsque le cardinal $k$ de $S$ est très grand. L'algorithme de Propp-Wilson propose une solution qui demande de simuler à chaque étape $k$ trajectoires de la chaîne. L'algorithme semble donc à priori inutilisable... Cependant, nous allons voir que si la chaîne considérée possède certaines propriétés de monotonie alors on peut réduire considérablement le nombres de trajectoires à simuler (en appliquant le principe du sandwich). L'algorithme devient alorsintéréssant. Nous allons commencer par un cas instructtif (mais sans aucun intérêt pratique) puis passer à un cas réellement intéressant (modèle d'Ising).\\

\dexer On désire simuler la probabilité uniforme $\pi$ sur $\{1,\ldots,k\}$ avec l'algorithme de Propp-Wilson en appliquant le principe du sandwich.\\
\quest 1) Montrer que la marche aléatoire simple symétrique sur $\{1, \ldots, k\}$ avec barrière réfléchissante en $1$ et en $k$ admet $\pi$ pour probabilité stationnaire.\\
\quest 2) Monter que la fonction $\phi$ définie par:
$$\phi(1,\cdot)=\indic_{[0,1/2)} + 2.\indic_{[1/2, 1]}$$
$$\phi(k,\cdot)=(k-1).\indic_{[0,1/2)} + k.\indic_{[1/2, 1]}$$
$$\phi(i,\cdot)=(i-1).\indic_{[0,1/2)} + (i+1).\indic_{[1/2, 1]}, \qquad i\in\{2,\ldots, k-1\}$$
est une fonction de mise à jour.\\
\quest 3) Montrer que la fonction $\phi$ vérifie
$$i \leq j \Rightarrow \phi (i, \cdot) \leq \phi (j, \cdot)$$
\quest 4) En déduire que toute trajectoire démarrant dans un état $i _in \{2, \ldots, k-1\}$ reste entre les trajectoires partant respectivement de $1$ et de $k$. Expliquer le terme ``principe du sandwich''.\\
\quest 5) Proposer une nouvelle version de l'algorithme de Propp-Wilson qui économise le nombre d'opérations.\\
\quest 6) Implémenter l'algorithme.
\fexer

\section{Le modèle d'Ising}

Soit $S$ un ensemble fini de $\ZZ^2$. Les point $s\in S$ sont appelés les sites. On introduit aussi l'espace des états $\Lambda=\{\pm 1\}$ et l'espace des configurations $E=\Lambda^S$.

\begin{deff}
Un champ aléatoire sur $S$ est une collection $X=\{X(s)\}_{s\in S}$ de variables aléatoires à valeurs dans $\Lambda$, où si l'on préfère une variable aléatoire à valeurs dans $\Lambda$.
\end{deff}

On dira qu'un site $y\in S$ est un voisin de $x\in S$ si et seulement si la distance (euclidienne standard) de $x$ à $y$ est égale à $1$ : $|x-y|=1$. On notera alors $x\sim y$. L'ensemble des voisins de $x$ dans $S$ est noté ${\cN}_x (S)=\{y \in S;\qquad |y-x|=1\}$ et on pose ${\tilde{\cN}}_x (S)= {\cN}_x (S) \cup \{x\}$. On remarquera que si $S=\ZZ^2$ alors pour tout $x$, le cardinal de ${\cN}_x (\ZZ^2)$ est $4$.

\begin{deff}
On dira que le champ aléatoire $X$ est un champ aléatoire Markovien sur $S$ si il vérifie : pour toute configuration  $x \in E$ et pour tout site $s\in S$
\begin{equation}
\PP (X(s)=x(s) | X(S\backslash s)= x(S \backslash s))= \PP (X(s)=x(s) | X({\cN}_s)=x({\cN}_s))
\end{equation}
Dans ce cas, la fonction $\pi^s: E \rightarrow [0,1]$ définie par $\pi^{s} (x)=\PP (X(s)=x(s) | X({\cN}_s)=x({\cN}_s))$ est appelée spécification locale du champ.
\end{deff}

Sous une condition de positivité assez faible sur la spécification locale, on peut montrer qu'il existe au plus une seule distribution de champ aléatoire Markovien pour une spécification locale donnée.

\subsection{Distribution de Gibbs}

\begin{deff}
On appelle clique sur $S$ tout singleton $\{s\}$ avec $s \in S$ ou bien toute paire $\{s,t\}$ de voisins $s,t \in S$. Une clique à deux éléments est dite maximale (ou bien non triviale).
\end{deff}

\begin{deff}
Un potentiel $V$ sur $E=\Lambda^S$ est une collection $\{V_C\}_{C\subset S}$ de fonctions de $E$ dans $\RR$ telles que :\\
i) $V_C= 0$ si $C$ n'est pas une clique.\\
ii) $V_C$ ne dépend de $x$ qu'à travers les valeurs de $x(C)$ (i.e. $V_C$ est à portée dans $C$).\\

La fonction énergie $\cE$ associé à un tel potentiel est définie par:
$$\cE (x)= \sum_{C \subset S} V_C (x), \qquad x \in E$$
Soit $T>0$. La distribution $\pi_T$ définie par
$$\pi_T (x)= Z_T^{-1} \exp\left(- \cfrac{\cE (x)}{T}\right), \qquad x \in E$$
où $Z_T$ est une constante de normalisation (dite fonction de partition) est appelée distribution de Gibbs à température $T$ associée au potentiel $V$.
\end{deff}

\noindent
\textit{Exemple : Modèle d'Ising.} C'est le modèle le plus simple que l'on puisse concevoir. On prend $S=\{1,\ldots, m\}^2$ et pour potentiel $V$,
$$V_{\{s\}} (x) = -Hx(s),\qquad V_{\{s,t\}}(x) = -J x(s) x(t) $$
où $H\in \RR$ est le champ magnétique externe et $J\in \RR^{+}$ est l'énergie interne d'un dipôle magnétique.

\begin{thf}
\label{thf-mes-gibbs}
Les distributions de Gibbs associées à un potentiel $V=\{V_{C}\}_{c \subset S}$ et à une température $T$ sont des champs Markoviens. Plus exactement, soit $X$ un champ aléatoire dont la loi est donnée par la distribution $\pi_T$ de la définition précédente, alors $X$ est un champ Markovien et ses spécifications locales sont données par :
\begin{equation}
\pi^s (x) =\cfrac{\exp (-T^{-1} \sum_{C \ni x} V_C (x))}{\sum_{ \lambda \in \Lambda} \exp (-T^{-1} \sum_{C \ni x} V_C (\lambda, x(S\backslash\{s\})}
\end{equation}
\end{thf}

\begin{remf}
La réciproque est quasi-vraie (il faut rajouter la condition de positivité dont on a parlé succintement). On pourra consulter $\cite{Br}$ pour une démonstration (où l'on verra une application inattendue de la formule de Möbius...).
\end{remf}

\subsection{Distribution de Gibbs en volume infini et transition de phase}

On cherche à étendre la notion de mesure de Gibbs lorsque $S$ est infini : $S=\ZZ^2$). Il est naturel, compte-tenu de la définition (\ref{thf-mes-gibbs}) de dire qu'une probabilité $\mu$ sur l'espace des configurations $\{\pm 1\}^{S}$ est une mesure de Gibbs associée au potentiel $V$ à température $T>0$ si elle vérifie les équations (dites DLR) suivantes :
$$\mu (X(s)=x(s) | X(S\backslash s)= x(S \backslash s))=\cfrac{\exp (-T^{-1}\sum_{C \ni x} V_C (x))}{\sum_{ \lambda \in \Lambda} \exp (-T^{-1}\sum_{C \ni x} V_C (\lambda, x(S\backslash\{s\})}$$
On dénote par $\cG (T)$ l'ensemble des mesures de Gibbs associée à ce potentiel à température $T$. Pour simplifier, on prendra le potentiel associé au modèle d'Ising avec le champ magnétique externe $H$ égal à $0$.\\
On cherche maintenant à caractériser les éléments de $\cG$. Pour cela, nous avons besoin d'introduire, en volume fini ($S$ fini), les mesures de Gibbs avec conditions aux bords qui généralisent les mesures de Gibbs introduites précedemment. Soit $S$ une ensemble fini et $W$ un sous-ensemble de $S$ (en général l'interieur de $S$). On se donne une configuration $\zeta$ sur $S\backslash W$, c'est-à-dire un élément de $\{0,1\}^{S\backslash W}$. La mesure de Gibbs $\pi_{T,\zeta}$ avec conditions aux bords $\zeta$ sur $S\backslash W$ à température $T>0$ est définie par :
 $$\forall \eta \in \{-1,1\}^{W}, \qquad \pi_{T,\zeta} (\eta)= Z_{T,\zeta}^{-1} \exp \left( -\cfrac{\cE (\eta^{\zeta})}{T}\right)$$
where $Z_{T,\zeta}$ est une fonction de normalisation et $\eta^{\zeta}$ est la configuration sur $\{-1,+1\}^{S}$ obtenue en concaténant $\eta$ et $\zeta$.\\
Deux cas particulièrement importants de mesures de Gibbs avec conditions aux bords peuvent être considérés. Prenons $S_m =\{1,\ldots,m\}^{2}$ où $m\geq 2$ et $W_{m}=\{2,\ldots,m-1\}^2$. On prend $\zeta$ la configuration sur $\{-1,1\}^{S_m\backslash W_m}$ composée uniquement de $+1$ (resp. de $-1$) et on note la mesure $\pi_{T,\zeta}$ obtenue $\pi^{+}_{T,m}$ (resp. $\pi^{-}_{T,m}$).

On peut alors démontrer le théorème suivant qui sert à caractériser l'ensemble des mesures de Gibbs en volume infini.
\begin{thf}
On considère le modèle d'Ising à température $T$.\\
1.  La limite (faible) de $(\pi^{+}_{T,m})_m$ (resp. de $\pi^{-}_{T,m})_m$) existe. On la note $\pi^+_T$ (resp. $\pi_T^{-}$).\\
2. $\pi_T^{\pm}$ sont des éléments de $\cG$.
\end{thf}

On dit qu'il y a transition de phases lorsque l'ensemble ${\cG}_T$ contient plus d'un élément. D'après le théorème précédent, on voit que si $\pi_T^{+} \neq \pi_T^{-}$, il y a transition de phase et on peut même démontrer la réciproque. On peut aller encore plus loin avec le théorème suivant.

\begin{thf}
Il y a transition de phase si et seulement si $\pi_T^{+} (\eta (0)=1) \neq \pi_T^{-} (\eta (0)=1)$. D'autre part il existe une température critique $T_c$ telle que si $T>T_c$, il n'y apas transition de phase et si $T<T_c$, il y a transition de phase.
\end{thf}

Finissons par quelques remarques. Nous avons ici considéré le cadre de la dimension $d=2$. Tout ceci se généralise en dimension $d$ quelconque. On peut démontrer que pour $d=1$, $T_c =0$ (il n'y a jamais transition de phase). Pour $d=2$, on connaît la valeur exacte de $T_c$, $T_c=2\left[\log(1+\sqrt{2})\right]$ (pour $J=1$). En dimension $d\geq 3$, c'est un problème ouvert. Le fait qu'à haute température, il y ait transition de phase est assez facile à montrer (argument de Peierls).




%\subsection{Transition de phase}
%
%On considère le modèle d'Ising cité en exemple dans la précédente section. $S$ est énuméré sous la forme $S= \{1, \ldots, N\}$ et avec les notations déjà introduites, on pose
%$$\pi_T (x) = Z_T^{-1} \exp \left( -T^{-1} {\cE} (x)\right)$$
%où
%$${\cE} (x) ={\cE}_0 (x) -H\sum_{i=1}^{N} x(i)$$
%avec
%$${\cE}_0 (x) = -\cfrac{J}{2} \sum_{s\sim t \in S} x(s) x(t) ={\cE}_0 (x)$$
%Le facteur $1/2$ apparu ici provient du fait que l'on compte maintenant les couples $(s,t)$ et non les paires $\{s,t\}$. Le terme $m(x)= \sum_{i=1}^{N} x( i)$ est appelé le moment magnétique de la configuration $x$. On introduit $M(H,T)=\pi_T (m)$ qui est le moment magnétique moyen lorsque le système est à l'équilibre sous la distribution de Gibbs $\pi_T$.
%
%\quest 1) Montrer que $-1 \leq M (H,T) \leq 1$.\\
%\quest 2) Montrer que $M(-H, T)= -M(H,T)$.\\
%\quest 3) Montrer que $\cfrac{\partial M}{\partial H} (H,T) \geq 0$.\\
%\quest 4) Montrer que $M (0,T)=0$.\\
%
%\subsection{Argument de Peierl

\subsection{Simulation d'une distribution de Gibbs par la méthode MCMC}

On considère un processus $(X_n)_{n \in \NN}$ à valeurs dans l'espace des configurations $E$. Autremant dit, pour tout $n \in \NN$, $X_n$ est une champ aléatoire.Supposons qu'il existe une chaîne de Markov irréductible aperiodique $(X_n)_{n \in \NN}$ à valeurs dans $E$ telle que sa distribution stationnaire soit donnée par$\pi_T$. Alors on a
\begin{equation}
\lim_{n\rightarrow \infty} d_V (\PP_{X_n}, \pi_T)=0
\end{equation}
où $d_V (\mu, \nu)$ est la distance en variation entre les deux probabilités $\mu$ et $\nu$ définie par
$$d_V (\mu , \nu)=\sum_{i \in E} |\mu(i) -\nu(i)|$$

On trouvera une preuve de ceci dans \cite{Br} ou \cite{Hag}. Cette propriété donne une idée de comment simuler une variable aléatoire $X$ ayant la probabilité $\pi_T$. La difficulté qui se présente lorsque l'on désire similer une variable aléatoire ayant pour loi $\pi_T$ est que l'on ne connaît pas la fonction de partition $Z_T$ et que le calcul de celle-ci, lorsque $E$ est très grand (de l'ordre de $2^{|S|}$), s'avère coûteux en temps et en précision.

L'échantillonage de Gibbs est le suivant. C'est en fait une version de l'algorithme de Metropolis. On se donne une probabilité $(q_s)_{s \in S}$ strictement positive (par exemple $q_s= 1/|S|$,i.e. la distribution un iforme sur $S$). La transition de probabilité qui permet de passer de $X_n =x$ à $X_{n+1}=y$ est obtenu de la manière suivante. On choisit un site $s \in S$ avec probabilité $q_s$. $x$ est alors changée en $y$ uniquement sur le site $s$, i.e. $x (S\backslash\{s\})=y(S\backslash \{s\})$; sur le site $s$ $x(s)$ est changé en $y(s)=\pm 1$ avec la probabilité
$$\pi (y(s)| x(S\backslash s))$$

La matrice de transition de la chaîne de Markov $(X_n)_n$ ainsi construite est
$$\PP (X_{n+1} =y | X_n =x) =q_s \pi (y (s) | x(S \backslash s)) \indic_{y(S\backslash s)= x(s \backslash s)}$$

La chaîne correspondante est irréductible et aperiodique (pourquoi ?). Pour cette chaîne, $\pi_T$ est la distribution stationnaire car elle satisfait la condition de ``bilan détaillé''. Le problème beaucoup plus compliqué est de savoir à quel moment (i.e. pour quel $n$) on peut considérer que la loi de $X_n$ est qusi égale à $\pi_T$. On consultera encore $\cite{BrG}$ à ce sujet. Un échantillonage plus astucieux existe. Il consiste, au lieu de ne changer que la valeur d'un site à chaque unité de temps, de changer la valeur des $|S|$ sites (cf. problème 7.6.2 de \cite{Br}).

\subsection{Simulation via l'algorithme de Propp-Wilson pour le modèle d'Ising}
On considère le modèle d'Ising à température $\beta^{-1}>0$. L'espace des configurations est $\Omega_m=\{-1,+1\}^{\Lambda_m}$ où $\Lambda_m=\{1,\ldots,m\}^2$. La probabilité que l'on considère sur l'espace des configurations est donnée par:
$$\mu_{m,\beta}(\zeta)=Z_{m,\beta}^{-1}\exp\left(\beta \sum_{x\sim y}\zeta(x) \zeta (y)\right)$$
où $Z_{m,\beta}$ est une constante de normalisation (\textit{fonction de partition}). Si $x$ est un site de $\Lambda_m$, on dénote par $k_{+} (x,\zeta)$ (resp. $k_{-} (x,\zeta)$) le nombre de voisins de $x$ qui ont pour valeur $+1$ (resp. $-1$).

\quest 1) Montrer que si $\zeta_{0}$ est un élément de $\Omega_m$ fixé
$$\mu_{m,\beta} (\zeta (x)=1 | \zeta(\Omega_m\backslash \{x\})=\zeta_0(\Omega_m \backslash \{x\}))=\cfrac{\exp\left(2\beta(k_{+}(x,\zeta_0)-k_{-}(x,\zeta_0)\right)}{\exp\left(2\beta(k_{+}(x,\zeta_0)-k_{-}(x,\zeta_0)\right)+1}$$

\quest 2) On introduit un ordre partiel sur $\Omega_m$ donné par $\zeta_0 \leq \eta_0$ si et seulement si pour tout site $x$, $\zeta_0 (x) \leq \eta_0 (x)$. Montrer que:
$$\mu_{m,\beta} (\zeta (x)=1 | \zeta(\Omega_m\backslash \{x\})=\zeta_0(\Omega_m \backslash \{x\}) ) \leq \mu_{m,\beta} (\zeta (x)=1 | \zeta(\Omega_m\backslash \{x\})=\eta_0(\Omega_m \backslash \{x\}))$$

\quest 3) Montrer qu'il existe une configuration maximale et une configuration minimale pour cet ordre.

\quest 4) Soit $(U_n)_{n}$ une suite de variables aléatoires uniformes sur $[0,1]$ indépendantes. On définit le processus de markov $(X_n)$ sur $\Omega_m$ de la manière suivante. Pour passer de $X_n$ à $X_{n+1}$, on choisit un site $x$ de $\Lambda_m$ de manière aléatoire et uniforme. On laisse inchangés les sites de $X_n$ autres que $x$. $X_{n+1} (x)$ est égal à $1$ si
$$U_{n+1} <\cfrac{\exp\left(2\beta(k_{+}(x,X_n)-k_{-}(x,X_n)\right)}{\exp\left(2\beta(k_{+}(x,X_n)-k_{-}(x,X_n)\right)+1}$$
et $-1$ sinon.

\subquest a) Montrer que $(X_n)_n$ admet $\mu_{m,\beta}$ pour probabilité stationnaire et que la chaîne est aperiodique et irréductible.

\subquest b) Monter que le principe du sandwich dans l'algorithme de Propp-Wilson s'applique pour cette chaîne.

\subquest c) Simuler une variable aléatoire de loi $\mu_{m,\beta}$ via l'algorithme de Propp-Wilson.

\subsection{Simulation de la transition de phase}

On propose d'illustrer par simulation informatique la transition de phase pour le modèle d'Ising $H=0$ et $J=1$.

\quest 1) Soit $m \geq 1$ et $T>0$. A l'aide de l'algoritme de Propp-Wilson, écrire une procédure qui simule une variable aléatoire ayant pour loi $\pi_{T,m}^{+}$ (resp. $\pi_{T,m}^{-}$) sur le carré $\{1,\ldots, m\}^2$. On note cette fonction $sim^{+} (m,T)$ (resp $\pi^{-} (m,T)$). Ce sont des tableaux de longueur $m^2$.

\quest 2) Ecrire une fonction $mom-magn^+ (T,m)$ (resp. $mom-magn^{-} (T,m)$) qui approxime $\pi_{T,m}^{+} (\zeta (0))$ (resp. $\pi_{T,m}^{-}(\zeta (0))$). Pour calculer cette espérance, on fera appel un certain nombre de fois à $sim^{+} (m,T)$ (resp. $sim^{-} (m,T)$) et on effectuera une moyenne empirique. $mom-magn^{\pm} (T,m)$ aura pour sortie un nombre.

\quest 3) Sur un graphique, tracer $mom-magn^{\pm} (m,T)$ en fonction de $T$. On choisira $m$ de taille assez grande pour considérer que $m\rightarrow \infty$ et pas trop grand pour que les ismulations ne durent pas plus de 10 minutes.

\quest 4) Observer que les deux courbes sont approximativement égales pour $T>T_c$ à la valeur $1/2$ mais pas pour $T<T_c$.

\quest 5) Attention, les simulations peuvent être assez longues. En particulier, on constatera que l'algorithme prend beaucoup plus de temps lorsque $T$ est aux environs de $T_c$. On évitera donc de prendre des valeurs trop proches de $T_c$.


\section{Méthode de gradient stochastique}

Une application à des méthodes de simulation apparaît dans ce que l'on appelle la méthode de gradient stochastique. Soit $E$ un ensemble fini et $U : E \rightarrow \RR$ une fonction appelée fonction de coût et que l'on cherche à minimiser. Autrement dit, on cherche à déterminer les $i_0$ de $E$ tels que
$$U(i_0)\leq U(i) \text{  pour tout $i\in E$}$$
Un tel élément est appelé un minimum global de la fonction de coût. L'algorithme de de gradient classique (ou déterministe) est le suivant et est grosso modo, une version discrète de la méthode du gradient. Pour chaque état $i\in E$, on se donne un ensemble (petit) $N(i)$ de sites de $E$ appelés les voisins de $i$. Supposons qu'à une certaine étape, une solution $i$ est examinée. A l'étape suivante, un site $j$ est choisi dans l'ensemble $N(i)$ des voisins de $i$ et on compare $U(i)$ avec $U(j)$. Si $U(i)$ est plus petit, on s'arrête là. Si $U(j)$ est plus petit, alors on repart de nouveau avec $j$ qui joue maintenant le rôle de $i$. Puisque $E$ est finie, l'algorithme s'arrêtera. Cependant, rien n'exclut qu'en fait de minimum global, ce soit un minimum local que l'on obtienne (c'est-à-dire que l'on tombe sur $i$ tel que $U(i) \leq U(j)$ pour tout $j\in N(i)$).\\
L'idée de l'algorithme de gradient stochastique est de donner la possibilité (avec une petite probabilité) à chaque étape de ne pas choisir le minimum entre $U(i)$ et $U(j)$. On évite ainsi de rester coincé sur un minimum local.\\

On se donne une suite de voisinages $N(i)$ qui communiquent :\\
Pour tout couple $i,j\in E$ d'états, il existe une trajectoire de $i$ à $j$, c'est-à-dire une suite finie d'états $i_1,\ldots,i_m \in E$ tels que $i_1\in N(i), i_2 \in N(i_1),\ldots, j\in N(i_m)$.\\

On définit maintenant une matrice de transition irréductible $(Q_{i,j})_{i,j \in E}$ telle que $Q_{i,j}>0$ si et seulement si $i=j$ ou $j\in N(i)$. On supposera de plus que $Q$ est symmétrique.\\

L'échantillonage de Métropolis est le suivant. On se donne un paramètre $T>0$ et on pose
$$\alpha_{i,j} (T)=\inf\left\{1,e^{\cfrac{(U(i)-U(j)}{T}}\right\}$$
On définit d'autre part une matrice de transition irréductible $P(T)$ par
$$P_{i,j}(T)=Q_{i,j}\alpha_{i,j}(T)$$

\quest 1) Montrer que si $U$ n'est pas une constante alors $P(T)$ est irréductible et aperiodique.

\quest 2) Montrer que $P(T)$ admet une probabilité stationnaire donnée par :
$$\forall i \in E, \; \; \pi_i (T)=\cfrac{e^{-U(i)/T}}{\sum_{k \in E}e^{-U(k)/T}}$$

\quest 3) Soit $H$ l'ensemble des minimas globaux
$$H=\{i \in E;\; \forall j\in E,\; U(i)\leq U(j)\}$$
Montrer que
\begin{equation*}
\lim_{T\rightarrow 0} \pi_{i}(T)=
\begin{cases}
\cfrac{1}{|H|} \text{  si $i\in H$}\\
0 \text{  si $i\notin H$}
\end{cases}
\end{equation*}

Ces résultats suggèrent la procédure suivante. On initialise l'algorithme avec une valeur du paramètre $T=T_0$ et on simule la chaîne de Markov de matrice de transition $P(T_0)$. Au bout d'un temps suffisament long, la distribution de la chaîne sera proche de la probabilité $\pi (T_0)$. A l'étape suivante, on choisit une valeur du paramètre $T=T_1 < T_0$ et à nouveau on simule la chaîne de Markov de matrice de transition $P(T_1)$ partant de la distribution obtenue (qui est proche de $\pi(T_0)$ donc proche de $\lim_{T\rightarrow 0} \pi(T)$). On attend un temps suffisament long pour que la distribution de la chaîne simulée soit proche de la distribution stationnaire $\pi (T_1)$. On continue ainsi de suite avec une suite de paramètres $(T_k)_{k\geq 0}$ qui tend vers $0$.\\
Toute la difficulté consiste alors à savoir quel temps $a_k$ il faut attendre entre deux changements de paramètres $T_k\rightarrow T_{k+1}$ pour être proche de la distribution stationnaire et quel doit être l'ordre de grandeur du passage de $T_k$ à $T_{k+1}$. Une condition nécessaire et suffisante a été donnée par Hajek (1988) mais elle est en pratique inutilisable. Elle énonce qu'il existe une constante $\gamma$ telle que l'algorithme précédent converge quelque soit l'état initial si et seulement si
$$\sum_{k=1}^{\infty} e^{-\cfrac{\gamma}{a_k}} =\infty$$
En particulier, on peut choisir $a_k$ de la forme
$$a_k=\cfrac{a}{\log(k+1)}$$
pour $a\geq \gamma$. Bien entendu, on ne connaît pas la valeur de $\gamma$.\\

Un exemple élémentaire (Exemple 13.4 p. 104-105 de \cite{Hag}) montre que si l'on choisit un refroidissement (i.e. une décroissance des $T_k$) trop rapide, alors l'algorithme peut defaillir. On considère un ensemble $E$ à $4$ éléments $s_1,\ldots, s_4$ et une fonction $f$ définie par:
$$f(s_1)=1, \qquad f(s_2)=2, \qquad f(s_3)=0,\qquad f(s_4)=2$$
On veut appliquer la méthode du gradient stochastique pour trouver le minimu de $f$.

\quest 1) Montrer que la matrice de transition de la chaîne de Markov qui permet d'appliquer l'algorithme pour une température $T$ est donnée par:
\begin{equation}
\left(
\begin{array}{cccc}
1-e^{-1/T} & 0.5 e^{-1/T}&0&0.5 e^{-1/T}\\
0.5&0&0.5&0\\
0&0.5e^{-2/T}&1-e^{-2/T}&0.5 e^{-2/T}\\
0.5&0&0.5&0
\end{array}
\right)
\end{equation}

Soit $(T_k)$ la suite des températures que l'on choisit dans l'algorithme et $a_k$ les temps de simulation avec température $T_k$. On commence donc par prendre une valeur $X_0$. On simule une chaîne de Markov partant de $X_0$ avec la température $T_0$ dans l'algorithme de Métropolis et au bout d'un certain temps $a_0$, on arrête là avec la température $T_0$. On a un nouvel état $X_1$. On part de ce nouvel état et on simule via Metropolis une chaîne de Markov mais à la température $T_1$ pendant $a_1$ et opn continue ainsi de suite. On suppose que l'on choisit $s_1$ pour valeur de $X_0$. Soit $A$ l'événement :''pour tout $k$, $X_k=s_1$.

\quest 1) Montrer que
$$\PP(A)=\prod_{i=1}^{\infty} (1-e^{-1/T_{i}})^{a_i}$$

\quest 2) Supposons que l'on choisisse d'appliquer l'algorithme avec $a_i = t_{0}=10^{5}$. Montrer que si $a_{i}=i^{-1}$ alors l'algorithme ne fonctionnera pas toujours.\\

\textit{Remarque :} Dans \cite{Hag}, il semble que $a_i =1$, ce qui est très étrange... Néanmoins, cela ne modifie en rien l'argumentation finale.
