3.2. Algorithmes de Metropolis

Soit \(\mu\) une mesure positive sur \(E\), non nulle telle que \(0 < \mu(x) < +\infty\) pour tout \(x \in E\). Les algorithmes type Metropolis permettent de construire une chaîne de Markov de transition \(Q\), vérifiant la condition de Doeblin, et dont la mesure \(\mu\) est invariante pour \(Q\), i.e. \(\mu Q = \mu\). Ainsi pour \(N \ge 1\) grand, la loi de \(X_N\) sera proche de la probabilité invariante \(\pi\) (la mesure \(\mu\) renormalisée, \(\pi(x) = \mu(x) / \sum_{y} \mu(y)\), et on pourra écrire un estimateur de Monte Carlo de la forme

\[\begin{equation*} \lim_{n \to +\infty} \frac{1}{n} \sum_{k=1}^n f(X_N^{(k)}) = \mathbf{E} \big[ f(X_N) \big] \simeq \sum_{x \in E} f(x) \pi(x), \end{equation*}\]

où \((X_N^{(k)})_{k \ge 1}\) est une suite i.i.d. de copies de \(X_N\). Un tel estimateur est appelé MCMC, et il y a 2 paramètres \(n \ge 1\) qui contrôle l’erreur statistique et \(N \ge 1\) qui peut être vu comme un paramètre de biais qui contrôle l’erreur commise en remplaçant \(\pi\) par la loi de \(X_N\).

3.2.1. Algorithme de Hastings Metropolis

Soit \(P\) une matrice stochastique irréductible telle que

(3.6)\[ \forall x, y \in E, \quad P(x, y) > 0 \quad \Leftrightarrow \quad P(y, x) > 0.\]

La matrice \(P\) est appelée matrice de proposition (ou matrice de sélection) et elle est utilisée pour parcourir l’espace d’état \(E\). On dit que deux éléments \(x\) et \(y\) sont voisins si \(P(x, y) > 0\).

Soit \(\rho: E \times E \to [0,1]\) la fonction définie par

(3.7)\[ \forall x, y \in E, \quad \rho(x, y) = \frac{\mu(y) P(y, x)}{\mu(x) P(x, y)} \wedge 1 = \min\Big\{\frac{\mu(y) P(y, x)}{\mu(x) P(x, y)}, 1 \Big\}.\]

La fonction \(\rho\) est le taux d’acceptation du voisin \(y\) proposé (en partant de la position \(x\)).

Voici l’algorithme de Mestropolis-Hastings que l’on va étudier dans la suite.

Algorithm 3.1 (Metropolis-Hastings)

Inputs Une mesure \(\mu\) (pas forcément normalisée), une matrice de proposition \(P\), un état initial \(x_0 \in E\) et un horizon \(N \ge 1\)

Output Une réalisation \((x_0, x_1, \dots, x_N)\) d’une chaine de Markov de matrice de transition \(Q\) invariante par \(\mu\)

• on définit \(\rho\) la fonction d’acceptation définie en (3.7)

• pour \(n \in \{0, \dots, N-1\}\)

  1. on propose un voisin \(y\) de \(x_n\) avec probabilité \(P(x_n, y)\) donnée par la matrice de proposition

  2. avec probabilité \(\rho(x_n, y)\) on accepte le voisin \(y\), et on pose \(x_{n+1} = y\), sinon on refuse l’état voisin et on pose \(x_{n+1} = x_n\).

Cet algorithme peut s’écrire sous la forme récursive suivante. Soit \((Y_n^{(x)})_{n \ge 1, x \in E}\) une suite de variable aléatoires indépendantes telle que

\[\begin{equation*} \forall x, y \in E, \quad \mathbf{P}\big(Y_n^{(x)} = y\big) = P(x, y) \end{equation*}\]

et une suite indépendante \((V_n)_{n \ge 1}\) i.i.d. telle que \(V_n \sim \mathcal{U}([0,1])\). La suite \((Y_n^{(x)})_{n \ge, x \in E}\) est doublement indicée et \(Y_n^{(x)}\) est le voisin proposé à l’itération \(n\) (pour construire \(x_n\)) lorsqu’on est dans l’état \(x \in E\) en \(n-1\) (i.e. \(x_{n-1} = x\)). Soit \(X_0\) indépendante des deux suites \(Y\) et \(V\), on définit

(3.8)\[ X_{n+1} = Y_{n+1}^{(X_n)} \mathbf{1}_{V_{n+1} \le \rho(X_n, Y_{n+1}^{(X_n)})} + X_n \mathbf{1}_{V_{n+1} > \rho(X_n, Y_{n+1}^{(X_n)})}.\]

Proposition 3.3

Le processus \((X_n)_{n \ge 0}\) défini par (3.8) est une chaine de Markov de matrice de transition \(Q\) définie par

\[\begin{equation*} Q(x, y) = \begin{cases} P(x, y) \rho(x, y) & \text{si $x \neq y$} \\ 1 - \sum_{z \neq x} Q(x, z) & \text{sinon} \end{cases} \end{equation*}\]

De plus, \((X_n)_{n \ge 0}\) est irréductible, réversible par rapport à \(\mu\) et donc \(\mu\) est invariante pour \(Q\).

Proof. D’après l’écriture récursive (3.8), on a \(X_{n+1} = \varphi(X_n, \xi_{n+1})\) où \(\varphi(x, z)\) est une fonction déterministe et \((\xi_n)_{n \ge 1}\) est une suite de variables aléatoires indépendantes, et indépendante de \(X_0\). La suite \((\xi_n)_{n \ge 1}\) est l’innovation donnée par \(\xi_n = (Y_n^{(x)}, x \in E, V_n)\). Donc \((X_n)_{n \ge 0}\) est une CM de noyau

\[\begin{align*} \forall x, y \in E, \quad Q(x, y) &= \mathbf{P}\big(X_{n+1} = y| X_n = x\big), \\ & =\mathbf{P}\big(Y_{n+1}^{(x)} \mathbf{1}_{V_{n+1} \le \rho(x, Y_{n+1}^{(x)})} + x \mathbf{1}_{V_{n+1} > \rho(x, Y_{n+1}^{(x)})} = y \big), \\ & = \begin{cases} \mathbf{P}\big(Y_{n+1}^{(x)} \mathbf{1}_{V_{n+1} \le \rho(x, Y_{n+1}^{(x)})} = y \big) & \text{si $x \neq y$,}\\ 1 - \sum_{z \neq x} Q(x, z) & \text{sinon.} \end{cases} \end{align*}\]

Or par indépendance on a

\[\begin{equation*} \mathbf{P}\big(Y_{n+1}^{(x)} \mathbf{1}_{V_{n+1} \le \rho(x, Y_{n+1}^{(x)})} = y\big) = \mathbf{P}\big(Y_{n+1}^{(x)} = y, V_{n+1} \le \rho(x, y) \big) = P(x, y) \rho(x, y). \end{equation*}\]

De plus, comme \(\mu(y) > 0\) pour tout \(y \in E\) et que \(P(y, x) > 0\) si et seulement si \(P(x, y) > 0\) on a \(\rho(x, y) > 0\) et l’irréductibilité de \(P\) implique l’irréductibilité de \(Q\). En effet s’il existe un chemin de longueur \(n\) reliant \(x\) et \(y\) i.e. un chemin \(x, x_1, \dots, x_n, y\) (d’éléments tous distincts) tel que \(P(x, x_1) \dots P(x_n, y) > 0\), alors

\[\begin{equation*} Q^{n}(x, y) \ge Q(x, x_1) \dots Q(x_n, y) = \rho(x, x_1) P(x, x_1) \dots \rho(x_n, y) P(x_n, y) > 0. \end{equation*}\]

Reste à montrer que \(Q\) est réversible par rapport à \(\mu\). Pour tout \(x, y \in E\),

\[\begin{equation*} \mu(x) Q(x, y) = \mu(x) P(x, y) \rho(x, y) = \mu(x) P(x, y) \wedge \mu(y) P(y, x). \end{equation*}\]

La quantité de droite est symétrique en \(x\) et \(y\) donc on a

\[\begin{equation*} \mu(x) Q(x, y) = \mu(y) Q(y, x). \end{equation*}\]

La proposition suivante est importante en pratique car elle permet de justifier que la loi de \(X_N\) est proche de \(\pi\) la probabilité invariante pour un \(N\) bien choisi (\(\pi\) est la mesure invariante \(\mu\) renormalisée \(\pi(x) = \mu(x) / \sum{y \in E} \mu(y)\)).

Proposition 3.4

On suppose \(E\) fini. Si \(P\) vérifie la condition de Doeblin, alors \(Q\) la vérifie aussi et donc \((X_n)_{n \ge 0}\) converge en variation vers \(\pi\) son unique probabilité invariante.

Proof. Soit

\[\begin{equation*} a = \min_{x, y \in E} \frac{\mu(y) P(y, x)}{\mu(x) P(x, y)} \mathbf{1}_{P(x,y) > 0} \le 1. \end{equation*}\]

Comme \(E\) est fini, on a \(a > 0\). D’autre part d’après la condition (3.6) on a par symétrie \(a \le 1\).
Comme \(\rho(x, y) \ge a\) on a

\[\begin{equation*} \forall x, y \in E, \, x \neq y, \quad Q(x, y) \ge a P(x, y). \end{equation*}\]

Montrons que cette minoration est vraie pour \(x = y\) i.e. \(Q(x, x) \ge a P(x, x)\). On a

\[\begin{equation*} Q(x, x) = 1 - \sum_{z \neq x} Q(x, z) \ge 1 - \sum_{z \neq x} P(x, z) = P(x, x), \end{equation*}\]

et donc \(Q(x, x) \ge a P(x, x)\) (car \(a \le 1\)). Ainsi

\[\begin{equation*} \forall x, y \in E, \quad Q(x, y) \ge a P(x, y), \end{equation*}\]

et pour tout \(n \ge 1\), \(Q^n(x, y) \ge a^n P^n(x, y)\). Si \(P\) vérifie la condition de Doeblin avec \(l \ge 1\), \(\alpha > 0\) et une probabilité \(c\) sur \(E\), alors \(Q\) vérifie la condition de Doeblin avec le même \(l\), la même probabilité \(c\) et \(\tilde \alpha = a^l \alpha\).

3.2.2. Algorithme de Metropolis (forme symétrique)

C’est la version originale de l’algorithme. Il s’agit de considérer une matrice de sélection \(P\) symétrique c’est à dire qu’on suppose \(P\) irréductible et

\[\begin{equation*} \forall x, y \in E, \quad P(x, y) = P(y, x). \end{equation*}\]

Cette condition est plus forte que la condition (3.6), et dans ce cas la fonction \(\rho(x, y)\) se simplifie en \(\rho(x, y) = \frac{\mu(y)}{\mu(x)} \wedge 1\) et la chaine de Markov \((X_n)_{n \ge 0}\) donnée par l’algorithme de Metropolis est de matrice de transition \(Q\) définie par

()\[\begin{equation} Q(x, y) = \begin{cases} \Big( \frac{\mu(y)}{\mu(x)} \wedge 1 \Big) P(x, y) & \text{si $x \neq y$}, \\ 1 - \sum_{z \neq x} Q(x, z) & \text{si $x = y$}. \end{cases} \end{equation}\]

On remarque immédiatement que si \(\mu\) est constante i.e. si \(\pi\) est la probabilité uniforme sur \(E\), alors \(\mu\) est réversible donc invariante pour \(P\) et \(Q = P\). On suppose donc que \(\mu\) n’est pas constante sur \(E\).

Proposition 3.5

On suppose \(E\) fini et \(\mu(x) > 0\), \(\forall x \in E\), non constante. Alors \(Q\) est irréductible et apériodique. La CM de matrice \(Q\) vérifie donc la condition de Doeblin et converge en variation vers \(\pi\).

Proof. Il suffit de prouver l’apériodicité, et pour cela il suffit de montrer qu’il existe \(x \in E\) tel que \(Q(x, x) > 0\). Intuitivement, si on propose un état \(y\) et qu’on le refuse alors on reste en \(x\) et donc \(Q(x, x) > 0\). On refuse un état avec probabilité \(1 - \rho(x, y)\) donc s’il existe \((x,y)\) tels que \(\rho(x, y) < 1\) i.e. \(\mu(y) < \mu(x)\) le résultat est prouvé. Comme la mesure \(\mu\) est non constante le résultat est vrai.

On peut aussi faire le raisonnement suivant: supposons que pour tout \(x \in E\), \(Q(x, x) = 0\). Alors

\[\begin{equation*} Q(x, x) - P(x, x) = \sum_{z \neq x} \big(1 - \rho(z, x)\big) P(x, z). \end{equation*}\]

La quantité à gauche est négative et celle de droite est positive donc les deux quantités sont nulles. On réécrit cette équation en échangeant \(z\) et \(x\) et on obtient

\[\begin{equation*} \sum_{z \neq x} \Big( \big(1 - \rho(z, x)\big) + \big(1 - \rho(x, z) \big) \Big) P(x, z) = 0. \end{equation*}\]

La quantité de gauche est strictement positive donc c’est absurde.