explicites autour du modèle d'Ising, on pourra consulter le livre [Bax82] de Rodney Baxter. Des modèles inspirés de celui d'Ising sont utilisés en imagerie pour la modélisation du bruit spatial.

La partie sur les algorithmes de Metropolis-Hastings, du recuit simulé, et de Propp-Wilson est inspirée d'un livre précédent [BC07, Chapitre 4] et du cours de Thierry Bodineau à l'École Polytechnique [Bod14]. L'algorithme de Metropolis-Hastings a été introduit par Nicholas Metropolis, Arianna Rosenbluth, Marshall Rosenbluth, Augusta Teller, et Edward Teller [MRR+53] puis généralisé par Keith Hastings [Has70], et fait aujourd'hui partie des méthodes MCMC (Monte Carlo Markov Chains), qui consistent à utiliser des chaînes de Markov pour approcher des espérances par la méthode de Monte Carlo. Les méthodes MCMC sont au cœur des implémentations quantitatives de la statistique bayésienne abordée par exemple dans le livre de Christian Robert et George Casella [RC04]. Des versions à particules permettent d'améliorer les performances pratiques notamment lorsque la loi $\mu$ est multimodale.

En métallurgie, le procédé du recuit consiste à recuire le métal pour échapper aux minima locaux d'énergie et obtenir une structure métallique de basse énergie, garantissant une meilleure solidité. L'algorithme du recuit simulé s'en inspire, ce qui explique son nom. Le recuit simulé converge théoriquement en temps infini lorsque la température suit un schéma de décroissance par paliers logarithmiques, mais ce résultat d'analyse asymptotique n'est pas vraiment pertinent en pratique. Une analyse de l'algorithme se trouve par exemple dans le livre de Étienne Pardoux [Par07] et dans le panorama de Olivier Catoxi [Cat99]. Des versions améliorées du recuit simulé, comme l'algorithme de Wang-Landau introduit par Fugao Wang et David Landau [WL01], constituent un sujet de recherche actuel. Les algorithmes de Metropolis-Hastings, du recuit simulé, et de Wang-Landau sont tous disponibles dans des cadres à temps et espace continus. L'algorithme de Propp-Wilson, proposé par James Propp et David Wilson [PW96, PW98], est présenté dans le livre de David Levin, Yuval Peres, et Elizabeth Wilmer [LPW09], et dans celui de Olle Häggström [Häg02]. L'exemple étudié dans la section 5.5 est inspiré d'un texte de Christophe Sabot. D'autres algorithmes de simulation exacte de la loi invariante ont été développés, comme par exemple celui de James Fill [Fil98].
<figure><img src="image_20.png" /></figure>
# Agrégation limitée par diffusion interne

**Mots-clés.** Modèle de croissance; ensembles aléatoires.

**Outils.** Chaîne de Markov; martingale; couplage; inégalité de concentration.

## Difficulté. *

Un fût de déchets radioactifs est enterré secrètement. Après quelques années, il devient poreux et laisse échapper son contenu. Pour éviter une contamination excessive, on a disposé des pièges à particules qui capturent la première particule qui passe mais deviennent ensuite inertes (une nouvelle particule passe sans être arrêtée). Le milieu est isotrope : une particule radioactive se déplace de la même manière dans toutes les directions tant qu'elle passe sur des pièges qui ont déjà été activés et est capturée par le premier piège libre qu'elle rencontre. On souhaite connaître la forme typique des zones qui seront contaminées par cette fuite. On introduit pour ce faire un modèle de croissance aléatoire, appelé *agrégation limitée par diffusion interne*, dont on étudie notamment les propriétés asymptotiques.

## 6.1 Modèle d'agrégation par diffusion

Pour modéliser l'évolution de la région polluée, on introduit le modèle de croissance suivant. On assimile l'espace à $Z^d$ où $d$ vaut 1, 2 ou 3 selon le nombre de dimensions que l'on souhaite prendre en compte. Deux éléments $x$ et $y$ de $Z^d$ sont dits voisins, et nous noterons $x \sim y$, si

$$|x - y|_1 = |x_1 - y_1| + \cdots + |x_d - y_d| = 1.$$

On place le fût à l'origine $0 \in Z^d$. Sur chaque point du réseau est disposé un piège. Notons $A(0)$ le singleton $\{0\}$. Le premier piège activé sera l'un des $2d$
emplacements voisins de l'origine. On modélise la trajectoire d'une particule radioactive sortant du fût par une marche aléatoire simple symétrique sur $Z^d$ issue de 0 et arrêtée lorsqu'elle quitte le point $A(0)$. On note $A(1)$ l'ensemble aléatoire composé de 0 et du point où la particule est sortie. Pour tout $x$ voisin de 0, $A(1)$ est égal à $\{0, x\}$ avec probabilité $1/2d$. On itère ensuite le procédé. Étant donné un ensemble $A(n) \subset Z^d$, on considère une marche aléatoire symétrique $(S_k)_{k \ge 0}$ issue de 0 arrêtée lorsqu'elle sort de $A(n)$. On définit alors $A(n+1)$ comme l'ensemble des éléments de $A(n)$ et du point où est sortie la marche $S$. Ce modèle est connu sous le nom de *agrégation limitée par diffusion interne*<fnref n="1" />.

On obtient ainsi une suite $(A(n))_{n \ge 0}$ d'ensembles aléatoires. Le résultat suivant regroupe quelques propriétés immédiates de cette suite.

**Théorème 6.1 (Propriétés immédiates).** La suite $(A(n))_{n \ge 0}$ est croissante au sens de l'inclusion. Le cardinal de $A(n)$ vaut exactement $n+1$. Pour tout $n \in \mathbb{N}$, soit $x$ et $y$ deux éléments de $A(n)$, alors il existe $x_0, \dots, x_m$ éléments de $A(n)$ tels que $x_0 = x$, $x_m = y$ et $x_i \sim x_{i+1}$ pour $i = 0, \dots, m-1$. On pourrait dire que l'ensemble $A(n)$ est connexe par arcs dans $Z^d$.

On souhaite à présent étudier le comportement asymptotique de la suite $(A(n))_{n \ge 0}$. On étudie en détail le cas de la dimension 1. L'étude en dimension supérieure, beaucoup plus délicate, est évoquée dans les sections 6.6 et 6.7.

## 6.2 Modèle unidimensionnel

On suppose que la propagation se fait selon un axe horizontal. On se place donc dans le cas où $d=1$. L'ensemble initial $A_0$ est le singleton $\{0\}$ puis $A_1$ est égal à $\{0, 1\}$ avec probabilité $1/2$ et $\{-1, 0\}$ avec probabilité $1/2$, etc. Notons $G_n = \min A_n$ et $D_n = \max A_n$. L'ensemble $A_n$ est de la forme $A_n = \{G_n, G_n+1, \dots, D_n-1, D_n\}$. Puisque le cardinal de $A_n$ est $n+1$, on a $D_n - G_n = n$. Ainsi, $A_n$ est caractérisé par $X_n = D_n + G_n$ et, en particulier,

$$D_n = \frac{X_n + n}{2} \quad \text{et} \quad G_n = \frac{X_n - n}{2}.$$

Les accroissements $(X_{n+1} - X_n)_{n \ge 0}$ ne peuvent prendre que les valeurs $-1$ ou $1$. Plus précisément, $X_{n+1} = X_n - 1$ si la marche issue de 0 atteint $G_n - 1$ avant $D_n + 1$, et $X_{n+1} = X_n + 1$ si la marche issue de 0 atteint $D_n + 1$ avant $G_n - 1$. La symétrie de la marche aléatoire sous-jacente fait que la loi de $X_n$ est symétrique pour tout $n \ge 0$. En particulier, pour tout $n \ge 0$, on a $0 = E(X_n) = E(D_n) + E(G_n)$ et donc $E(D_n) = n/2$ et $E(G_n) = -n/2$.

**Théorème 6.2 (Évolution markovienne).** La suite $(X_n)_{n \ge 0}$ est une chaîne de Markov inhomogène sur $Z$ issue de 0 dont les transitions sont décrites par les relations suivantes : pour $-n \le i \le n$,

<footnote n="1">1. « Internal Diffusion Limited Aggregation (IDLA) » en anglais.</footnote>
$$P(X_{n+1} = i - 1 | X_n = i) = \frac{n+2+i}{2(n+2)}$$

et

$$P(X_{n+1} = i + 1 | X_n = i) = \frac{n+2-i}{2(n+2)}.$$

*Démonstration.* Pour tout $n \ge 1$,

$$\{X_n = i\} = \left\{ D_n = \frac{i+n}{2}, G_n = \frac{i-n}{2} \right\}.$$

La probabilité $P(X_{n+1} = i - 1 | X_n = i)$ est donc égale à la probabilité que la marche aléatoire simple $(S_n)_{n \ge 0}$ issue de 0 atteigne $(i-n)/2 - 1$ avant $(i+n)/2 + 1$. Or d'après le théorème 2.2 (ruine du joueur) si $a \le 0$ et $b \ge 0$ sont deux entiers distincts, et si $T_i := \inf \{n \ge 0, S_n = i\}$ pour $i \in \{a, b\}$ et $T := T_a \wedge T_b$, alors, $T$ est intégrable et

$$P(T = T_a) = 1 - P(T = T_b) = \frac{b}{b-a} \quad \text{et} \quad E(T) = -ab.$$

□

**Remarque 6.3 (Équilibrage).** Lorsque $X_n$ est strictement positif, $X_{n+1} - X_n$ a une probabilité plus forte de valoir -1 que de valoir +1. La tendance s'inverse lorsque $X_n$ est strictement négatif. On peut donc penser que le processus aura davantage tendance à revenir en 0 qu'une marche aléatoire simple symétrique (dont l'accroissement vaut plus ou moins 1 avec probabilité 1/2).

## 6.3 Convergence presque sûre

On établit ici le premier résultat de convergence en temps long pour le processus $(X_n)_{n \ge 0}$ du modèle unidimensionnel. L'idée de la preuve est de rendre rigoureuse l'intuition de la remarque 6.3 en comparant la valeur absolue de $(X_n)_{n \ge 0}$ à celle de la marche aléatoire simple par une méthode de couplage.

**Théorème 6.4 (Convergence presque sûre).** Le processus $(|X_n|)_{n \ge 0}$ est une chaîne de Markov inhomogène et

$$\frac{|X_n|}{n} \xrightarrow{n \to \infty} 0 \quad \text{d'où} \quad \frac{G_n}{n} \xrightarrow{n \to \infty} -\frac{1}{2} \quad \text{et} \quad \frac{D_n}{n} \xrightarrow{n \to \infty} \frac{1}{2}.$$

*Démonstration.* En discutant suivant les valeurs possibles de $X_n$, on constate que $|X_{n+1}| = 1$ si $|X_n| = 0$, et que si $|X_n| > 0$ alors

$$|X_{n+1}| = \begin{cases} |X_n| - 1 & \text{avec probabilité} \quad \frac{1}{2} + \frac{|X_n|}{2(n+2)}, \\ |X_n| + 1 & \text{avec probabilité} \quad \frac{1}{2} - \frac{|X_n|}{2(n+2)}. \end{cases}$$
Ainsi $|(X_n)|_{n \ge 0}$ est une chaîne de Markov (suite récurrente aléatoire) inhomo-
gène. Pour établir que $|X_n|/n \to 0$ presque sûrement, on procède par couplage.
Plus précisément, on construit deux processus $(Y_n)_{n \ge 0}$ et $(Z_n)_{n \ge 0}$ comme
suit : $Y_0 = Z_0 = 0$ et pour $n \ge 0$,

$$
Y_{n+1} = \begin{cases} Y_n + 2\mathbf{1}_{\{U_{n+1} < (n+2-Y_n)/2(n+2)\}} - 1 & \text{si } Y_n > 0, \\ 1 & \text{si } Y_n = 0, \end{cases}
$$

et

$$
Z_{n+1} = \begin{cases} Z_n + 2\mathbf{1}_{\{U_{n+1} < 1/2\}} - 1 & \text{si } Z_n > 0, \\ 1 & \text{si } Z_n = 0, \end{cases}
$$

où $(U_n)_{n \ge 1}$ est une suite de v.a.r. i.i.d. de loi uniforme sur $[0, 1]$. Le processus $(Y_n)_{n \ge 0}$ a la loi de $|(X_n)|_{n \ge 0}$ tandis que le processus $(Z_n)_{n \ge 0}$ a même loi que la valeur absolue d'une marche aléatoire simple. Ce couplage est tel que presque sûrement, pour tout $n \in \mathbb{N}$, $0 \le Y_n \le Z_n$. En effet, si $0 < Y_n \le Z_n$ alors, par construction, $Y_{n+1} \le Z_{n+1}$. Par ailleurs, si $0 = Y_n \le Z_n$ alors, $Z_n$ peut être nul et dans ce cas, $Y_{n+1} \le Z_{n+1}$ ou $Z_n \ge 2$ car $Z_n$ et $Y_n$ ont même parité. Enfin, $Z_n/n \to 0$ p.s. grâce à la loi forte des grands nombres. $\square$

**Remarque 6.5** (Espérance). Le théorème 6.2 donne, avec $a_{n+1} := \frac{n+1}{n+2}$,

$$
E(|X_{n+1}| | |X_n|) = \mathbf{1}_{|X_n|=0} + a_{n+1} |X_n| \mathbf{1}_{|X_n| \neq 0}
$$

d'où, en observant que $|X_n|\mathbf{1}_{|X_n|\ne 0} = |X_n|$,

$$
E(|X_{n+1}|) = \mathbb{P}(X_n = 0) + a_{n+1}E(|X_n|).
$$

6.4 Convergence en loi

On s'intéresse à présent aux fluctuations de $X_n/n$ autour de sa limite.
Dans cette section, on donne des éléments de preuve du résultat suivant.

**Théorème 6.6 (Convergence en loi).** La convergence en loi suivante a lieu :

$$
\frac{X_n}{\sqrt{n}} \xrightarrow{n \to \infty} \mathcal{N}\left(0, \frac{1}{3}\right).
$$

**Remarque 6.7 (Comparaison avec le théorème limite central).** La vitesse de convergence en loi de $(X_n)_{n \ge 0}$ est la même que pour la marche aléatoire simple et la loi limite est encore gaussienne. Pourtant sa variance asymptotique est plus petite que dans le cas identiquement distribué. Cela provient de la tendance auto-stabilisatrice de la chaîne $(X_n)_{n \ge 0}$.

Avant d'aborder la démonstration du théorème 6.6, étudions le compor-
tement des premiers moments de $X_n/n$ quand $n$ tend vers l'infini. La loi de
$X_n$ est symétrique donc tous ses moments impairs sont nuls. Pour le moment d'ordre 2, on a, pour $n \ge 0$,

$$
\begin{align*}
\mathrm{E}[X_{n+1}^2 | \mathcal{F}_n] &= \frac{n+2-X_n}{2(n+2)}(X_n+1)^2 + \frac{n+2+X_n}{2(n+2)}(X_n-1)^2 \\
&= 1 + \frac{n}{n+2}X_n^2.
\end{align*}
$$

En posant $x_n(2) = \mathrm{E}(X_n^2)$ et en prenant l'espérance dans la relation ci-dessus,
on obtient, pour $n \ge 1$,

$$
x_{n+1}(2) = 1 + \frac{n}{n+2}x_n(2) \quad \text{et} \quad x_1(2) = 1.
$$

Une récurrence donne

$$
x_n(2) = 1 + \frac{1}{n(n+1)} \sum_{k=1}^{n-1} k(k+1),
$$

et $(x_n(2)/n)_{n \ge 1}$ converge donc vers 1/3. On peut établir de même que pour tout $k \ge 0$, la suite $(x_n(2k)/n^k)_{n \ge 1}$ admet une limite $\mu(2k)$ et que la suite $(\mu(2k))_{k \in \mathbb{N}}$ est solution de

$$
\mu(0) = 1 \quad \text{et} \quad \mu(2k+2) = \frac{2k+1}{3}\mu(2k).
$$

Or d'après le théorème 21.7, la seule mesure de probabilité dont les moments
pairs vérifient la récurrence ci-dessus et les moments impairs sont nuls est
la loi $\mathcal{N}(0, 1/3)$. D'autre part, d'après le théorème 21.6, la convergence des
moments vers une mesure de probabilité caractérisée par ses moments entraîne
la convergence en loi. On obtient ainsi la convergence en loi de $X_n/\sqrt{n}$ vers
la mesure gaussienne centrée de variance 1/3.

Un autre moyen d'établir la convergence annoncée est d'utiliser un raisonnement basé sur le théorème limite central pour les martingales. Il existe plusieurs jeux d'hypothèses plus ou moins faciles à vérifier. Voici le théorème que nous utiliserons dans la suite.

**Théorème 6.8 (TLC pour martingales).** *Soit $(M_n)_{n \ge 0}$ une martingale de carré intégrable pour une filtration $(\mathcal{F}_n)_{n \ge 0}$, et soit $(a_n)_{n \ge 0}$ une suite réelle, strictement positive, déterministe croissante vers l'infini. On pose*

$$
\langle M \rangle_0 = 0 \quad \text{et} \quad \langle M \rangle_{n+1} = \langle M \rangle_n + \mathrm{E}((\Delta M_{n+1})^2 | \mathcal{F}_n),
$$

où $\Delta M_n := M_n - M_{n-1}$. Supposons que

1. il existe $\lambda \ge 0$ déterministe tel que $a_n^{-1} \langle M \rangle_n \xrightarrow{n \to \infty} \lambda$,

2. (condition de Lyapunov) il existe $\delta > 0$ tel que

$$
\frac{1}{(a_n)^{1+\delta/2}} \sum_{k=1}^{n} \mathrm{E}(\Delta M_k^{2+\delta} | \mathcal{F}_{k-1}) \xrightarrow{n \to \infty} 0.
$$
Alors, on a

$$ \frac{1}{\sqrt{a_n}} M_n \xrightarrow{n \to \infty} \mathcal{N}(0, \lambda), \quad \text{et} \quad \sqrt{a_n} \frac{M_n}{\langle M \rangle_n} \xrightarrow{n \to \infty} \mathcal{N}(0, \lambda^{-1}) \quad \text{si } \lambda > 0. $$

Il convient de modifier légèrement $(X_n)_{n \ge 0}$ pour en faire une martingale.

**Lemme 6.9** (Martingale dans le modèle). Le processus $(M_n)_{n \ge 0}$ défini par $M_n = (n+1)X_n$ est une martingale de carré intégrable pour la filtration naturelle $(\mathcal{F}_n)_{n \ge 0}$ de $(X_n)_{n \ge 0}$, qui vérifie

$$ \mathrm{E}((\Delta M_{n+1})^2 | \mathcal{F}_n) = (n+2)^2 - X_n^2, $$

et

$$ \mathrm{E}((\Delta M_{n+1})^4 | \mathcal{F}_n) = ((n+2)^2 - X_n^2)((n+2)^2 + 3X_n^2), $$

et

$$ \frac{\langle M \rangle_n}{n^3} \xrightarrow{n \to \infty} \frac{1}{3}. $$

**Démonstration.** Soit $\alpha \in \{1, 2, 4\}$. La quantité suivante

$$ \begin{aligned} \mathrm{E}((\Delta M_{n+1})^{\alpha} | \mathcal{F}_n) &= \mathrm{E}(((n+2)X_{n+1} - (n+1)X_n)^{\alpha} | X_n) \\ &= [X_n - (n+2)]^{\alpha} \frac{n+2+X_n}{2(n+2)} + [X_n + n + 2]^{\alpha} \frac{n+2-X_n}{2(n+2)}, \end{aligned} $$

est nulle si $\alpha = 1$, et donc $(M_n)_{n \ge 0}$ est une martingale. Si $\alpha \in \{2, 4\}$ alors

$$ \mathrm{E}((\Delta M_{n+1})^{\alpha} | \mathcal{F}_n) = ((n+2)^2 - X_n^2) \frac{(n+2-X_n)^{\alpha-1} + (n+2+X_n)^{\alpha-1}}{2(n+2)}, $$

ce qui donne les deux premières relations attendues. En particulier, on a

$$ \langle M \rangle_n = \sum_{k=0}^{n-1} ((k+2)^2 - X_k^2) = \frac{n^3}{3} - \sum_{k=0}^{n-1} X_k^2 + o(n^3). $$

De plus

$$ \frac{1}{n^3} \sum_{k=0}^{n-1} X_k^2 \le \frac{1}{n} \sum_{k=1}^{n-1} \frac{X_k^2}{k^2}. $$

Puisque $X_n/n \to 0$ p.s. le lemme de Cesàro permet de conclure. $\square$

<figure><img src="image_21.png" /></figure>

Pour pouvoir obtenir la convergence en loi de $(M_n)_{n \ge 0}$ correctement renormalisée, il reste à vérifier la condition de Lyapunov du théorème 6.8 avec $a_n = n^3$. Pour $\delta = 2$, on a

$$ \sum_{k=1}^{n} \mathrm{E}((\Delta M_k)^4 | \mathcal{F}_{k-1}) \le \sum_{k=2}^{n+1} k^4 + 2 \sum_{k=1}^{n} (k+1)^2 X_{k-1}^2 = \mathcal{O}(n^5). $$

La suite $(M_n/n^{3/2})_{n \ge 1}$ converge donc en loi vers $\mathcal{N}(0, 1/3)$. Ceci fournit immédiatement la convergence dans le théorème 6.6.
## 6.5 Une inégalité de concentration

Complétons les résultats sur le comportement de $X_n$ par une inégalité de concentration valable pour tout $n \ge 1$ (à la différence du théorème de convergence en loi qui n’est qu’asymptotique).

**Théorème 6.10 (Concentration).** Pour tout $n \ge 1$ et tout $r > 0$,

$$ P\left(\frac{|X_n|}{n} \ge r\right) \le 2 \exp\left(-\frac{3n^2r^2}{8(n+1)}\right). $$

*Démonstration.* Par définition de la martingale $(M_n)_{n \ge 0}$ définie dans le lemme 6.9, pour tout $r > 0$,

$$ P\left(\frac{|X_n|}{n} \ge r\right) \le P(|M_n| \ge n(n+1)r). $$

La martingale $(M_n)_{n \ge 0}$ a des accroissements borné car pour tout $n \ge 1$,

$$ |M_n - M_{n-1}| = |(n+1)X_n - nX_{n-1}| \le n|X_n - X_{n-1}| + |X_n| \le 2n $$

puisque $|X_n - X_{n-1}| = 1$ et $|X_n| \le n$. D'après l'inégalité de concentration d'Azuma-Hoeffding du lemme 20.3, pour tout $\lambda > 0$, on a

$$ P(|M_n| \ge \lambda) \le 2 \exp \left( - \frac{\lambda^2}{2 \sum_{k=1}^{n} (2k)^2} \right). $$

Grâce à la comparaison série-intégrale

$$ \sum_{k=1}^{n} (2k)^2 \le 4 \int_{1}^{n+1} x^2 dx = \frac{4}{3}(n+1)^3, $$

il reste à remplacer $\lambda$ par $n(n+1)r$ pour obtenir l'inégalité souhaitée. $\square$

## 6.6 Dimensions supérieures

On se place à présent sur $\mathbb{Z}^d$ avec $d > 1$. Comme en dimension 1, on s'attend à ce que la forme limite de l'ensemble aléatoire devienne relativement régulière puisque la marche aléatoire issue de l'origine aura tendance à atteindre l'ensemble aléatoire $A_n^c$ en des points qui ont beaucoup de voisins dans $A_n$. La figure 6.1 fournit une réalisation de $A(1000)$ en dimension 2.

Pour tout $r \ge 0$ notons $\|x\|$ la norme euclidienne d'un point $x \in \mathbb{R}^d$ et

$$ B(0, r) = \{y \in \mathbb{R}^d : \|y\| < r\} \text{ et } \mathbb{B}(0, r) = B(0, r) \cap \mathbb{Z}^d. $$

Le cardinal d'un ensemble $\Lambda$ de $\mathbb{Z}^d$ est noté $|\Lambda|$. Fixons un entier $n$ et considérons l'ensemble aléatoire $A(|\mathbb{B}(0, n)|)$. Le rayon (aléatoire) de la plus grande boule incluse dans cet ensemble est
<figure><img src="image_22.png" /><figcaption>Fig. 6.1. Une réalisation de A(1000) en dimension 2.</figcaption></figure>

$$ \sup \{r \ge 0 : \mathbb{B}(0, r) \subset A(|\mathbb{B}(0, n)|)\} $$

tandis que celui de la plus petite boule qui le contient est

$$ \inf \{r \ge 0 : A(|\mathbb{B}(0, n)|) \subset \mathbb{B}(0, r)\} $$

On définit alors l'erreur interne $\delta_I(n)$ comme

$$ n - \delta_I(n) = \sup \{r \ge 0 : \mathbb{B}(0, r) \subset A(|\mathbb{B}(0, n)|)\}, $$

et l'erreur externe $\delta_E(n)$ comme

$$ n + \delta_E(n) = \inf \{r \ge 0 : A(|\mathbb{B}(0, n)|) \subset \mathbb{B}(0, r)\}. $$

Le théorème suivant montre non seulement que la forme asymptotique de $A(n)$ est sphérique mais fournit également une borne supérieure pour les fluctuations extrêmes de $A(n)$.

**Théorème 6.11** (Forme et fluctuation). La forme limite de $(A_n)_{n \ge 0}$ est la trace sur $\mathbb{Z}^d$ d'une boule euclidienne. De plus, il existe une constante $\beta_d$ telle que presque sûrement

— pour $d \ge 3$ :

$$ \limsup_{n \to \infty} \frac{\delta_I(n)}{\sqrt{\log(n)}} \le \beta_d \quad \text{et} \quad \limsup_{n \to \infty} \frac{\delta_E(n)}{\sqrt{\log(n)}} \le \beta_d, $$
— pour $d = 2$ :

$$ \limsup_{n \to \infty} \frac{\delta_I(n)}{\log(n)} \le \beta_2 \quad \text{et} \quad \limsup_{n \to \infty} \frac{\delta_E(n)}{\log(n)} \le \beta_2. $$

## 6.7 Pour aller plus loin

Le modèle d'agrégation limitée par diffusion interne<fnref n="2" /> est considéré par Persi Diaconis et William Fulton dans [DF91]. Sur le versant de la modélisation, il peut être relié à des méthodes de polissage chimique : des particules sont lâchées dans une cuve et détruisent un élément du bord rendant la cuve plus lisse. On trouvera dans les articles [AG13a, AG13b] d'Amine Asselah et d'Alexandre Gaudillière des références dans cette direction.

L'une des première études de $(A(n))_{n \ge 0}$ est due à Maury Bramson, David Griffeath et Gregory Lawler dans [LBG92]. Il y est notamment établi que la forme asymptotique de $A(n)$ est sphérique. Il faudra attendre longtemps pour avoir une étude complète pour les fluctuations autour de cette forme limite. On pourra retrouver le théorème 6.11 dans les articles de David Jerison, Lionel Levin, et Scott Scheffied [JLS12, JLS13] et dans [AG13a, AG13b].

Le modèle en dimension 1 peut être relié à un modèle de boules et d'urnes dit de Bernard Friedman décrit dans l'article [Fre65] de David Freedman<fnref n="3" />. Une boule tirée est remise dans l'urne accompagnée de $\alpha$ boules de la même couleur et $\beta$ boules de la couleur opposée. Si $\beta = 0$, on retrouve le modèle de l'urne de Pólya étudié dans le chapitre 15. Lorsque $\beta > 0$, on peut montrer que la proportion asymptotique de chaque couleur tend vers 1/2 pour tout $\alpha$.

Il existe également un modèle plus difficile à étudier, connu sous le nom de *agrégation limitée par diffusion externe*, dans lequel les particules proviennent de l'infini et viennent s'agglutiner au fil du temps.

Le théorème limite central pour les martingales 6.8 se trouve dans le livre de Peter Hall et Christopher Heyde [HH80, Cor. 3.1], et dans [BC07, Th. 3.35].

<footnote n="2">2. «Internal Diffusion Limited Aggregation» (IDLA) en anglais.</footnote>

<footnote n="3">3. fried-man ≠ freed-man!</footnote>
<figure><img src="image_23.png" /></figure>
# Chaînes de Markov cachées

**Mots-clés.** estimation paramétrique; estimation récursive; variable cachée; espérance conditionnelle; processus autorégressif.

**Outils.** Chaîne de Markov; vecteur gaussien; martingales; maximum de vraisemblance; formule de Bayes.

## Difficulté. *

Bien que markoviens, de nombreux processus réels ne sont observés que partiellement. Ce chapitre présente deux modèles de *chaînes de Markov causées*, qui sont par définition des projections de chaînes de Markov. Le défi mathématique est alors d'obtenir à partir des observations partielles le plus d'informations possible sur la partie cachée. Il se ramène à un problème de calcul de lois conditionnelles. Il n'est en général pas possible de calculer explicitement ces lois. C'est cependant le cas dans les deux exemples présentés dans ce chapitre : l'un s'appuie sur les chaînes de Markov à espace d'états fini, l'autre sur un modèle gaussien.

## 7.1 Algorithme progressif-rétrograde

On modélise la structure d'un brin d'ADN 1 par une chaîne de caractères écrite dans un alphabet fini $\mathcal{A}$. Certaines parties de la chaîne correspondent à la suite des codes des acides aminés nécessaires à la fabrication d'une protéine, et on dit qu'elles sont codantes, tandis que d'autres parties de la chaîne ne le sont pas. Pour tenir compte de cette structure segmentée, on modélise la chaîne par une trajectoire de chaîne de Markov possédant deux régimes.

À un brin d'ADN de longueur $l$, noté $(X_1, \dots, X_l)$, et à valeurs dans $\mathcal{A}$, on adjoint le $l$-uplet $(U_1, \dots, U_l)$, à valeurs dans $\mathcal{U} = \{0,1\}$. Si la position

<footnote n="1">1. Acide DésoxyriboNucléique, vecteur du code génétique des organismes vivants.</footnote>
$n$ est codante (respectivement non codante) alors $U_n = 1$ (respectivement $U_n = 0$). On modélise la loi de $(X, U) = ((X_n, U_n))_{n \ge 0}$ par une chaîne de Markov homogène sur $\mathcal{A} \times \mathcal{U}$ de matrice de transition donnée par

$$ 
\begin{aligned}
P((x, u), (x', u')) &= \mathbb{P}(X_{n+1} = x', U_{n+1} = u' | X_n = x, U_n = u) \\
&= \rho(u, u')\pi_u(x, x'),
\end{aligned}
$$

où $\rho$ est une matrice de transition sur $\{0, 1\}$ et où $\pi_0$ et $\pi_1$ sont des matrice de transitions sur $\mathcal{A}$. On suppose que tous les coefficients des matrices $\rho$, $\pi_0$ et $\pi_1$ sont strictement positifs. La composante $U$ agit comme une sorte d'interrupteur ou de commutateur (markovien) à deux positions.

La chaîne $(X, U)$ est irréductible, récurrente et apériodique car tous les coefficients de sa matrice de transition sont strictement positifs. En général la projection $X = (X_n)_{n \ge 0}$, qui modélise le brin d'ADN, n'est pas une chaîne de Markov. Le modèle est rigide : la loi de la longueur d'un segment de la catégorie $u \in \mathcal{U}$ suit une loi géométrique de paramètre $1 - \rho(u, u)$.

On se place à présent dans la situation où la matrice de transition de $(X, U)$, c'est-à-dire le triplet $\rho, \pi_0, \pi_1$, est connue, par exemple grâce à une estimation menée au préalable (nous y reviendrons plus loin). La question est à présent la suivante : peut-on déterminer la loi de $(U_j)_{1 \le j \le l}$ conditionnellement aux observations $(X_i)_{1 \le i \le l}$ ? On dit que la chaîne de Markov $(X, U)$ est cachée car la composante $U$ n'est pas observée. Seule la composante $X$ est observée.

**Exemple 7.1.** Voici un exemple utilisé dans la figure 7.1 :

## Description de l'algorithmme

On adopte la notation vectorielle $X_{1:i}$ et $x_{1:i}$ pour désigner les $i$-uplets $(X_1, \dots, X_i)$ et $(x_1, \dots, x_i)$. On définit les probabilités suivantes, pour tous $1 \le i \le l$, $x \in \mathcal{A}^l$, et $v \in \mathcal{U}$ :

— Probabilité de prédiction :

$$ P^i(v) := \mathbb{P}(U_i = v | X_{1:i-1} = x_{1:i-1}). $$

C'est la probabilité que l'état caché $U_i$ soit en position $v$ connaissant les observations $X_{1:i-1}$;

— Probabilité de filtrage :

$$ F^i(v) := \mathbb{P}(U_i = v | X_{1:i} = x_{1:i}). $$

C'est la probabilité que l'état caché $U_i$ soit en position $v$ connaissant les observations $X_{1:i}$;
— Probabilité de lissage :

$$L^i(v) := \mathbb{P}(U_i = v | X_{1:l} = x_{1:l}).$$

C'est la probabilité que l'état caché $U_i$ soit en position $v$ connaissant les observations $X_{1:l}$.

Pour calculer les probabilités de lissage, on utilise un algorithme *progressif-rétrograde*<fnref n="2" /> : on détermine par récurrence, dans le sens croissant des indices, les probabilités de prédiction et de filtrage, puis on en déduit, toujours par récurrence mais descendante, les probabilités de lissage. Les relations de récurrence nécessaires à cet algorithme sont rassemblées dans le théorème suivant.

**Théorème 7.2** (Prédiction, filtrage, lissage). Pour tous $1 \le i \le l$ et $v \in \mathcal{U}$,

$$ 
\begin{align*}
P^i(v) &= \sum_{u \in \mathcal{U}} \rho(u, v) F^{i-1}(u), \\
F^i(v) &= \frac{\pi_v(x_{i-1}, x_i) P^i(v)}{\sum_{u \in \mathcal{U}} \pi_u(x_{i-1}, x_i) P^i(u)}, \\
L^{i-1}(u) &= F^{i-1}(u) \sum_{v \in \mathcal{U}} \rho(u, v) \frac{L^i(v)}{P^i(v)}.
\end{align*}
$$

*Démonstration.* Pour l'équation de prédiction, on fait intervenir $U_{i-1}$ puis d'utiliser le fait que les lois $\text{Loi}(U_i | U_{i-1}, X_{1:i-1})$ et $\text{Loi}(U_i | U_{i-1})$ coïncident :

$$ 
\begin{align*}
P^i(v) &= \sum_{u \in \mathcal{U}} \mathbb{P}(U_{i-1} = u, U_i = v | X_{i-} = x_{i-}) \\
&= \sum_{u \in \mathcal{U}} \mathbb{P}(U_i = v | U_{i-1} = u, X_{i-} = x_{i-}) \mathbb{P}(U_{i-1} = u | X_{i-} = x_{i-}) \\
&= \sum_{u \in \mathcal{U}} \rho(u, v) F^{i-1}(u)
\end{align*}
$$

avec la notation allégée $x_{i-} := x_{1:i-1}$. L'utilisation de la relation de base $\mathbb{P}(A|B \cap C)\mathbb{P}(B|C) = \mathbb{P}(A \cap B|C)$ permet de traiter l'équation de filtrage :

$$ 
\begin{align*}
F^i(v) &= \frac{\mathbb{P}(U_i = v, X_i = x_i | X_{1:i-1} = x_{i-})}{\mathbb{P}(X_i = x_i | X_{i-} = x_{i-})} \\
&= \frac{\mathbb{P}(U_i = v, X_i = x_i | X_{i-} = x_{i-})}{\sum_{u \in \mathcal{U}} \mathbb{P}(U_i = u, X_i = x_i | X_{i-} = x_{i-})} \\
&= \frac{\mathbb{P}(X_i = x_i | U_i = v, X_{i-} = x_{i-}) \mathbb{P}(U_i = v, | X_{i-} = x_{i-})}{\sum_{u \in \mathcal{U}} \mathbb{P}(X_i = x_i | U_i = u, X_{i-} = x_{i-}) \mathbb{P}(U_i = u | X_{i-} = x_{i-})}.
\end{align*}
$$

Or par définition de la matrice de transition de $(X, U)$,

$$\mathbb{P}(X_i = x_i | U_i = u, X_{1:i-1} = x_{1:i-1}) = \mathbb{P}(X_i = x_i | U_i = u, X_{i-1} = x_{i-1})$$

<footnote n="2">2. En anglais : « forward-backward ».</footnote>
$$= \pi_{u_{i-1}}(x_i - 1, x_i).$$

Pour l'équation de lissage, on écrit

$$ \mathbb{P}(U_{i-1} = u, U_i = v | X_{1:l} = x_{1:l}) = \mathbb{P}(U_{i-1} = u | U_i = v, X_{1:l} = x_{1:l}) \mathbb{P}(U_i = v | X_{1:l} = x_{1:l}), $$

puis

$$ \mathbb{P}(U_{i-1} = u | U_i = v, X_{1:l} = x_{1:l}) = \frac{\mathbb{P}(U_{i-1} = u, X_{1:i-1} = x_{1:i-1}, X_{i+1:l} = x_{i+1:l} | U_i = v, X_i = x_i)}{\mathbb{P}(X_{1:i-1} = x_{1:i-1}, X_{i+1:l} = x_{i+1:l} | U_i = v, X_i = x_i)}. $$

La propriété de Markov assure que, conditionnellement au présent, le passé et le futur sont indépendants. On obtient donc après simplification que

$$ \mathbb{P}(U_{i-1} = u | U_i = v, X_{1:l} = x_{1:l}) = \mathbb{P}(U_{i-1} = u | U_i = v, X_{1:i} = x_{1:i}). $$

On écrit, à nouveau, grâce à l'expression de la matrice de transition de $(X, U)$,

$$ \begin{aligned} \mathbb{P}(U_{i-1} = u | U_i = v, X_{1:i} = x_{1:i}) &= \frac{\mathbb{P}(U_{i-1} = u, U_i = v, X_i = x_i, | X_{1:i-1} = x_{1:i-1})}{\mathbb{P}(U_i = v, X_i = x_i | X_{1:i-1} = x_{1:i-1})} \\ &= \frac{F^{i-1}(u)\rho(u, v)\pi_v(x_{i-1}, x_i)}{P^i(v)\pi_v(x_{i-1}, x_i)}. \end{aligned} $$

On a donc obtenu

$$ \mathbb{P}(U_{i-1} = u, U_i = v | X_{1:l} = x_{1:l}) = F^{i-1}(u)\rho(u, v) \frac{L^i(v)}{P^i(v)}, $$

ce qui fournit la relation de lissage en sommant sur $v$. $\square$

On initialise l'algorithme en choisissant pour $P^1(u)$ la loi initiale (par exemple la loi stationnaire). Les relations de prédiction et filtrage du théorème 7.2 permettent de calculer toutes les probabilités $P^i$ et $F^i$ par récurrence progressive. On en déduit ensuite les probabilités $L^i$ par récurrence rétrograde grâce à la relation de lissage du théorème 7.2 avec l'initialisation $L^l(v) = F^l(v)$. La figure 7.1 illustre l'efficacité de l'algorithme : la plupart du temps la probabilité $L^i(v)$ est supérieure à 1/2 quand $U_i = v$ (dans 95,2% des cas exactement dans cet exemple). On voit toutefois dans cet exemple que l'algorithme rate une zone où $U$ vaut 1 (aux alentours de 120) et qu'il a toujours un peu de retard aux changements de régimes. L'algorithme fonctionne d'autant mieux que les matrices $\pi_0$ et $\pi_1$ sont différentes (les deux régimes ont des comportements très différents) et que $\rho(0, 0)$ et $\rho(1, 1)$ sont grands (les plages entre deux changements de régime sont longues).
<figure><img src="image_24.png" /><figcaption>Fig. 7.1. Trajectoire de U (courbe continue) et probabilité de lissage pour $v = 0$ (courbe en pointillés) avec $l = 200$ et les transitions de l'exemple 7.1.</figcaption></figure>

Estimation d'une matrice de transition

L'étude précédente suppose connues les matrices de transition. Il est toutefois important de savoir estimer la matrice de transition d'une chaîne de Markov à partir d'une de ses trajectoires. On considère une chaîne de Markov $(Z_n)_{n \ge 0}$ sur l'espace d'états fini $E$ de matrice de transition $\theta_* = (\theta_*(i, j))_{(i,j) \in E^2}$. On suppose que la chaîne est irréductible et récurrente, et on note $\mu_*$ sa loi invariante. On suppose la mesure initiale $\nu$ connue mais la matrice de transition $\theta_*$ inconnue. On note $\Theta$ l'ensemble des matrices de transition sur $E$. Soit $\lambda$ la mesure de comptage sur $E$. La loi de $(Z_0, \dots, Z_n)$ a pour densité par rapport $\lambda^{\otimes n+1}$ la fonction $(z_0, \dots, z_n) \in E^{n+1} \mapsto f_{\theta_*}(z_0, \dots, z_n)$ donnée par

$$f_{\theta_*(z_0, \dots, z_n)} = \mathbb{P}(Z_0 = z_0, \dots, Z_n = z_n) = \nu(z_0)\theta_*(z_0, z_1) \cdots \theta_*(z_{n-1}, z_n).$$

Pour tous $\theta \in \Theta$ et $z \in E^{n+1}$ on note

$$f_{\theta}(z_0, \dots, z_n) := \nu(z_0)\theta(z_0, z_1) \cdots \theta(z_{n-1}, z_n).$$

La quantité aléatoire $L_{Z_1, \dots, Z_n}(\theta) := f_{\theta}(Z_0, \dots, Z_n)$ est appelée vraisemblance<fnref n="3" /> de $\theta$ pour l'échantillon $(Z_0, \dots, Z_n)$. L'estimateur de maximum de vraisemblance $\hat{\theta}_n$ de $\theta_*$ est par définition l'élément aléatoire de $\Theta$ qui maximise la fonction aléatoire $\theta \in \Theta \mapsto L_{Z_0, \dots, Z_n}(\theta)$.

<footnote n="3">3. En anglais, on dit likelihood, d'où la notation $L$.</footnote>
**Lemme 7.3 (Maximum de vraisemblance).** *L'estimateur de maximum de vraisemblance $\hat{\theta}_n$ de $\theta_*$ est donné pour tous $i, j \in E$ par*

$$ \hat{\theta}_n(i,j) = \begin{cases} \frac{N_n^{ij}}{N_n^i} & \text{si } N_n^i > 0, \\ 0 & \text{sinon,} \end{cases} \quad \text{où} \quad \begin{cases} N_n^i := \sum_{p=0}^{n-1} 1_{\{Z_p=i\}}, \\ N_n^{ij} := \sum_{p=0}^{n-1} 1_{\{Z_p=i, Z_{p+1}=j\}}. \end{cases} $$

Notons que si $\theta(i, j) = 0$ alors $\hat{\theta}_n(i, j) = 0$, mais que la réciproque est fausse.

*Démonstration.* On numérote $E$ de sorte que $E = \{1, \dots, s\}$. Notons $L := L_{Z_1, \dots, Z_n}$. Il est plus commode de maximiser $\log(L)$ plutôt que $L$. On doit maximiser $\theta \in \Theta \mapsto \log(L(\theta))$ sous les contraintes $\theta(i, j) \ge 0$ et $\sum_{k=1}^s \theta(i, k) = 1$ pour tous $1 \le i, j \le s$. Or on a

$$ \partial_{\theta(i,j)}(\theta \mapsto \log(L(\theta))) = \frac{N_n^{ij}}{\theta(i,j)}, $$

et les $s$ conditions d'extrémalité reviennent à dire que pour tout $1 \le i \le s$ la fonction $j \mapsto N_n^{ij}/\theta(i, j)$ doit être constante. Ce qui donne le résultat en tenant compte de $\theta \in \Theta$. $\square$

L'estimateur $\hat{\theta}$ est intuitif puisqu'il estime $\theta(i, j)$ par la proportion de transitions de l'état $i$ vers l'état $j$ parmi les $n$ sauts observés. Les propriétés asymptotiques de cet estimateur se déduisent du comportement des suites $(N_n^{ij})_n$ et $(N_n^i)_n$ qui sont établies ci-dessous.

**Théorème 7.4 (Convergence et normalité asymptotique).** *Pour tout $x \in E$, sachant $\{X_0 = x\}$, l'estimateur $\hat{\theta}_n$ est convergent et asymptotiquement normal : pour tous $i, j \in E$,*

$$ \hat{\theta}_n(i, j) \xrightarrow{n \to \infty} \theta(i, j) $$

$$ \text{et} \quad \sqrt{n\mu(i)}(\hat{\theta}_n(i, j) - \theta(i, j)) \xrightarrow{n \to \infty} \mathcal{N}(0, \theta(i, j)(1 - \theta(i, j))). $$

*Démonstration.* On a, pour tous $x, i, j \in E$, sachant $\{X_0 = x\}$,

$$ \frac{1}{n} N_n^i \xrightarrow{n \to \infty} \mu(i) \quad \text{et} \quad \frac{1}{n} N_n^{ij} \xrightarrow{n \to \infty} \mu(i)\theta(i, j) $$

tandis que
$$ \frac{1}{\sqrt{n}}(N_n^{ij} - N_n^i \theta(i, j)) \xrightarrow{n \to \infty} \mathcal{N}(0, \mu(i)\theta(i, j)(1 - \theta(i, j))). $$

Ces convergences presque sûres découlent de la loi des grands nombres pour les chaînes de Markov $Z$ et $Y$, où $Y$ est la chaîne de Markov irréductible et récurrente d'espace d'états $E^2$ définie par $Y_n = (Z_n, Z_{n+1})$ et dont la loi invariante est $(i, j) \mapsto \mu(i)\theta(i, j)$. La convergence en loi découle du théorème limite central pour les chaînes de Markov. $\square$
## 7.2 Filtre de Kalman

Un avion se déplace entre Paris et Londres, en tentant de suivre une trajectoire théorique définie par le plan de vol. L'avion est surveillé au sol par des contrôleurs aériens grâce à un radar qui reçoit un écho de l'avion à intervalles réguliers. La trajectoire effective de l'avion s'écarte de la trajectoire théorique pour de multiples raisons (météorologie, imprécision du pilote automatique, turbulences,...). On cherche donc à localiser l'avion au cours de son vol à partir des observations radars successives.

On note $X_n$ l'écart (inconnu) entre la trajectoire théorique et la position de l'avion au temps $n$. De plus, on note $Y_n$ la mesure donnée par le radar au temps $n$. Cette mesure est entachée d'erreurs à cause de l'imprécision du radar. Le problème qui se pose à l'aiguilleur est d'estimer au mieux la position de l'avion au temps $n$ au vu des observations $Y_0, \dots, Y_n$. Pour simplifier l'étude, on supposera que l'objet observé évolue dans un espace de dimension 1. On modélise les suites aléatoires $(X_n)_{n \ge 0}$ et $(Y_n)_{n \ge 0}$ en posant

$$ egin{cases} X_0 = W_0, \\ X_n = aX_{n-1} + W_n \quad n \ge 1, \\ Y_n = X_n + V_n \quad n \ge 0, \end{cases} $$

où $a$ est un nombre réel déterministe, et où $(V_n)_{n \ge 0}$ et $(W_n)_{n \ge 0}$ sont des suites aléatoires indépendantes, avec $(V_n)_{n \ge 0}$ i.i.d. de loi gaussienne $N(0, \tau^2)$ et $(W_n)_{n \ge 0}$ i.i.d. de loi gaussienne $N(0, \sigma^2)$. Les variables aléatoires $(W_n)_{n \ge 0}$ représentent les fluctuations instantanées de l'écart entre position théorique et position réelle. Les variables aléatoires $(V_n)_{n \ge 0}$ modélisent les erreurs de mesure du radar. Le paramètre $a$ modélise l'action du pilote.

En théorie du signal, on dit que $(X_n)_{n \ge 0}$ est un processus autorégressif d'ordre 1, noté AR(1), à bruit gaussien, et on a

$$ X_n = \sum_{k=0}^{n} a^k W_{n-k} = \sum_{k=0}^{n} a^{n-k} W_k, \quad n \in \mathbb{N}. $$

À présent, le problème est d'estimer $X_n$ sachant les observations $Y_0, \dots, Y_n$. Le candidat naturel est l'espérance conditionnelle $E(X_n | Y_0, \dots, Y_n)$ puisque c'est la fonction des observations qui est la plus proche de $X_n$ dans $L^2$ (moindre carrés). Or le caractère gaussien du modèle permet, grâce aux propriétés des vecteurs gaussiens, de calculer explicitement la loi conditionnelle $\text{Loi}(X_n | Y_0, \dots, Y_n)$, en particulier sa moyenne qui est $E(X_n | Y_0, \dots, Y_n)$.

La première remarque importante est que pour tout $n \ge 1$,

$$ (X_0, \dots, X_n, Y_0, \dots, Y_n) $$

est un vecteur aléatoire gaussien puisque toute combinaison linéaire de ses coordonnées est une combinaison linéaire des variables aléatoires gaussiennes
indépendantes $(W_i)_{0 \le i \le n}$ et $(V_i)_{0 \le i \le n}$. En particulier, pour tout $n \ge 0$, $(X_n, Y_0, \dots, Y_n)$ est aussi un vecteur gaussien. De plus, la loi de $X_n$ sachant $(Y_0, \dots, Y_n)$ est une loi gaussienne, dont on note la moyenne $\hat{X}_n$ et la variance $P_n$. On se propose de calculer ces quantités par récurrence. Comme pour les chaînes de Markov à espace d'états fini de la première partie du chapitre, la récurrence se fait en deux étapes. L'étape de prédiction consiste à exprimer la loi $\text{Loi}(X_n | Y_0, \dots, Y_{n-1})$ en fonction de $\text{Loi}(X_{n-1} | Y_0, \dots, Y_{n-1})$. Puis, dans l'étape de filtrage, on prend en compte l'observation $Y_n$ pour exprimer $\text{Loi}(X_n | Y_0, \dots, Y_n)$ en fonction de $\text{Loi}(X_n | Y_0, \dots, Y_{n-1})$.

**Lemme 7.5 (Lois conditionnelles).** On a

$$ \text{Loi}(X_n | Y_0, \dots, Y_n) = \mathcal{N}(\hat{X}_n, P_n) $$

où

$$ \hat{X}_n = a\hat{X}_{n-1} + \frac{P_n}{\tau^2}(Y_n - a\hat{X}_{n-1}) \quad \text{et} \quad P_n = \frac{a^2\tau^2 P_{n-1} + \sigma^2\tau^2}{a^2P_{n-1} + \sigma^2 + \tau^2}. $$

*Démonstration.* Tout d'abord, on rappelle le résultat suivant sur les vecteurs gaussiens, que nous appelons formule de Bayes : si $(X, Y_0, \dots, Y_{n-1})$ est un vecteur gaussien dans $\mathbb{R}^{n+1}$ avec

$$ \text{Loi}(X | Y_0, \dots, Y_{n-1}) = \mathcal{N}(\mu, \gamma^2) \quad \text{et} \quad \text{Loi}(Y | Y_0, \dots, Y_{n-1}, X) = \mathcal{N}(X, \delta^2) $$

alors

$$ \text{Loi}(X | Y_0, \dots, Y_{n-1}, Y) = \mathcal{N}\left(\rho^2\left(\frac{\mu}{\gamma^2} + \frac{Y}{\delta^2}\right), \rho^2\right) \quad \text{où} \quad \frac{1}{\rho^2} := \frac{1}{\gamma^2} + \frac{1}{\delta^2}. $$

À présent, et comme annoncé, on procède par récurrence, en deux étapes.
— *Initialisation.* Puisque $Y_0 = X_0 + V_0$, on a $\text{Loi}(Y_0 | X_0) = \mathcal{N}(X_0, \tau^2)$.
La formule de Bayes assure que $\text{Loi}(X_0 | Y_0) = \mathcal{N}(\hat{X}_0, P_0)$ où

$$ \hat{X}_0 = \frac{\sigma^2}{\sigma^2 + \tau^2} Y_0 \quad \text{et} \quad P_0 = \frac{\sigma^2 \tau^2}{\sigma^2 + \tau^2}. $$

— *Prédiction.* On a $\text{Loi}(X_{n-1} | Y_0, \dots, Y_{n-1}) = \mathcal{N}(\hat{X}_{n-1}, P_{n-1})$, et grâce au modèle,

$$ \text{Loi}(X_n | Y_0, \dots, Y_{n-1}) = \mathcal{N}(a\hat{X}_{n-1}, a^2 P_{n-1} + \sigma^2). $$

— *Filtrage.* D'après le modèle à nouveau, il vient

$$ \text{Loi}(Y_n | Y_0, \dots, Y_{n-1}, X_n) = \mathcal{N}(X_n, \tau^2). $$

On applique alors la formule de Bayes pour inverser le conditionnement entre $Y_n$ et $X_n$: la loi $\text{Loi}(X_n | Y_0, \dots, Y_n) = \mathcal{N}(\hat{X}_n, P_n)$ avec

$$ \frac{1}{P_n} = \frac{1}{a^2 P_{n-1} + \sigma^2} + \frac{1}{\tau^2} \quad \text{et} \quad \hat{X}_n = P_n \left( \frac{a \hat{X}_{n-1}}{a^2 P_{n-1} + \sigma^2} + \frac{Y_n}{\tau^2} \right). $$
☐

**Remarque 7.6** (Gain suite aux observations). *La majoration*

$$P_n = \mathbb{E}[(X_n - \hat{X}_n)^2] \le \mathbb{E}[(X_n - Y_n)^2] = \tau^2,$$

peu surprenante connaissant les propriétés de l'espérance conditionnelle, sou-
ligne bien que l'on gagne effectivement à utiliser toutes les observations
$Y_0, \dots, Y_n$ plutôt que de se contenter de la dernière $Y_n$.

Estimation de certains paramètres du modèle

On suppose dans cette section qu'on observe les positions de l'avion sans erreurs, c'est-à-dire qu'on a accès à la suite $(X_i)_{1\le i\le n}$. On souhaite estimer les coefficients $a$ et $\sigma^2$. La question n'est pas complètement évidente car les observations ne sont pas indépendantes.

**Lemme 7.7** (Estimateur de maximum de vraisemblance). *L'estimateur de maximum de vraisemblance $(\hat{a}, \hat{\sigma})$ de $(a, \sigma)$ est donné par*

$$\hat{a}_n = \frac{\sum_{k=1}^{n} X_{k-1} X_k}{\sum_{k=1}^{n} X_{k-1}^2} \quad \text{et} \quad \hat{\sigma}_n^2 = \frac{1}{n} \sum_{k=1}^{n} (X_k - \hat{a}_n X_{k-1})^2.$$

*Démonstration*. Il suffit de maximiser le logarithme de la vraisemblance donnée par

$$L_{(X_1, \dots, X_n)}(a, \sigma) = \frac{1}{(2\pi\sigma^2)^{n/2}} \prod_{k=1}^{n} \exp \left( - \frac{(X_k - aX_{k-1})^2}{2\sigma^2} \right).$$

☐

**Théorème 7.8** (Convergence et normalité). *Si $-1 < a < 1$ alors*

$$(\hat{a}_n, \hat{\sigma}_n) \xrightarrow{n \to \infty} (a, \sigma),$$

et

$$\sqrt{n}(\hat{a} - a) \xrightarrow{n \to \infty} \mathcal{N}(0, 1 - a^2), \quad \text{et} \quad \sqrt{n}(\hat{\sigma}^2 - \sigma^2) \xrightarrow{n \to \infty} \mathcal{N}(0, 2\sigma^4).$$

*Démonstration*. On se contente d'établir les résultats sur $\hat{a}$, en supposant $\sigma$
connu et fixé. On réécrit $\hat{a}_n$ de la manière suivante :

$$\hat{a}_n = a + \frac{\sum_{k=1}^{n} X_{k-1} W_k}{\sum_{k=1}^{n} X_{k-1}^2} = a + \frac{M_n}{\sum_{k=1}^{n} X_{k-1}}$$

où $M_n = \sum_{k=1}^{n} X_{k-1}W_k$, avec $M_0 = 0$. La suite $(M_n)_{n \ge 0}$ est une martingale
par rapport à la filtration $(\mathcal{F}_n)_n$ naturelle de $W$, de processus croissant
$$ \langle M \rangle_0 = 0 \quad \text{et, pour } n \ge 1, \quad \langle M \rangle_n = \sigma^2 \sum_{k=1}^{n} X_{k-1}^2. $$

À présent, on vérifie tout d'abord, en utilisant le fait que $X$ est un processus autorégressif d'ordre 1 AR(1) avec $|a| < 1$, que

$$ \frac{\langle M \rangle_n}{n} \xrightarrow{n \to \infty} \frac{\sigma^4}{(1 - a^2)}, $$

(on a en particulier $\langle M \rangle_n \to +\infty$ p.s. quand $n \to \infty$), puis on utilise la loi des grands nombres et le théorème limite central pour les martingales de carré intégrable, qui donnent

$$ \frac{M_n}{\langle M \rangle_n} \xrightarrow{n \to \infty} 0 \quad \text{et} \quad \frac{M_n}{\sqrt{\langle M \rangle_n}} \xrightarrow{n \to \infty} \mathcal{N}(0, 1). $$

□

## 7.3 Pour aller plus loin

L'algorithme de segmentation progressif-rétrograde remonte au moins aux travaux des années 1960 de Leonard Baum et Lloyd Welch, et peut être vu comme une instance du concept général de programmation dynamique. Les chaînes de Markov cachées ont été utilisées notamment pour la reconnaissance de la parole dans les années 1970, et pour la génomique à partir des années 1980. Le livre de Stéphane Robin, François Rodolphe, et Sophie Schbath [RRS05] propose une introduction accessible à l'utilisation des chaînes de Markov cachées en génomique. On peut également consulter à ce sujet le livre de Étienne Pardoux [Par07]. Bien que les modèles utilisés en pratique soient plus sophistiqués que celui présenté dans ce chapitre, notamment en ce qui concerne les espaces d'état $\mathcal{A}$ et $\mathcal{U}$, ils font appel aux mêmes concepts et outils. Cependant en pratique, on ne connaît pas en général les matrices de transition ni même le nombre d'états cachés pertinent pour rendre compte de la loi de la séquence, et il faut donc construire des algorithmes qui permettent en plus d'estimer ces paramètres de complexité. On pourra également consulter le livre de Jean-François Delmas et Benjamin Jourdain [DJ06] à ce sujet. D'autre part, il est possible de tester si une suite aléatoire est markovienne ou pas en utilisant par exemple $N^i$ et $N^{ij}$ et le test du $\chi^2$, comme expliqué par exemple dans le livre de Didier Dacunha-Castelle et Marie Duflo [DCD83].

La partie sur le filtre de Kalman est inspirée du livre de David Williams [Wil91]. Le filtre de Kalman, présenté ici dans une version simple, est un grand classique de la théorie du signal, développé dès les années 1960 par Thorvald Thiele et Peter Swerling, et par Rudolf Kalman et Richard Bucy, notamment pour les besoins du programme Apollo de la National Aeronautics and Space Administration. On trouvera dans [Par07] l'expression du filtre de
Kalman en dimension supérieure à 1. La loi des grands nombres et le théorème limite central pour les martingales utilisé dans la preuve du théorème 7.8 sont tirés du livre [BC07]. L'estimation par maximum de vraisemblance des paramètres $a, \sigma, \tau$ du modèle gaussien à partir de $Y$ sans connaître $X$ se trouve par exemple dans le livre de Peter Brockwell et Richard Davis [BD02, Sec. 8.5].
<figure><img src="image_25.png" /></figure>
# Algorithme EM et mélanges

**Mots-clés.** Mélange de lois; estimation paramétrique.

**Outils.** Formule de Bayes; maximum de vraisemblance; entropie relative; loi gaussienne.

**Difficulté.** *

L'algorithme Expectation-Maximization (EM) fait partie des algorithmes les plus importants de la statistique. Il permet d'approcher numériquement l'estimateur de maximum de vraisemblance pour les modèles partiellement observés. De nombreuses variantes sont disponibles. Dans ce chapitre, nous illustrons la mécanique fondamentale de l'algorithme EM sur un exemple simple.

Sur une île cohabitent quatre espèces de mouettes différentes. Les ornithologues souhaitent estimer la proportion de mouettes de chaque espèce à partir de l'observation de la taille des nids. Contrairement aux oiseaux, les nids ne bougent pas, ce qui facilite le comptage. Malheureusement, les différentes espèces font des nids assez ressemblants : on ne sait pas par quelles espèces ils ont été construits. On suppose en revanche que la distribution des nids est caractéristique d'une espèce. La distribution globale de la taille des nids apparaît comme le mélange de quatre lois de probabilité, chacune rendant compte de la répartition de la taille des nids pour chaque espèce d'oiseaux.

## 8.1 Modèle pour les nids de mouettes

Modélisons la distribution de la taille des nids construits par l'espèce *j* par une loi $\mu_j = N(m(j), v(j))$ de densité $\gamma_{m(j),v(j)}$. La loi $\mu$ de la taille des nids de mouettes admet donc pour densité la fonction

$$x \mapsto \sum_{j=1}^{J} \alpha(j) \gamma_{m(j),v(j)}(x) = \sum_{j=1}^{J} \frac{\alpha(j)}{\sqrt{2\pi v(j)}} \exp \left( -\frac{(x - m(j))^2}{2v(j)} \right),$$
avec $\alpha_1 > 0, \dots, \alpha_J > 0$ et $\alpha_1 + \dots + \alpha_J = 1$. Cette combinaison convexe
finie de densités de probabilité est appelée mélange fini<fnref n="1" /> de populations. Si
$X$ désigne la taille du nid et $Z$ l'espèce de l'oiseau qui l'a construit, la loi du
couple $(X, Z)$ est donnée de la manière suivante :

— la variable Z suit la loi discrète sur {1, ..., J} de poids α(1), ..., α(J):

$$
\text{Loi}(Z) = \sum_{j=1}^{J} \alpha(j) \delta_j
$$

— pour tout $1 \le j \le J$, conditionnellement à $\{Z = j\}$, la variable $X$ suit
la loi $\mathcal{N}(m(j), v(j))$ :

$$
\mathrm{Loi}(X|Z=j) = \mathcal{N}(m(j), v(j)).
$$

Les paramètres $(\alpha_j, m(j), v(j))_{1\le j\le J}$ sont inconnus. Le vecteur $\alpha$ représente les coefficients du mélange tandis que les vecteurs $m$ et $v$ représentent respectivement les vecteurs des moyennes et des variances des lois gaussiennes. On souhaite estimer toutes ces quantités à la seule vue de la taille de $n$ nids. Pour cela, on suppose qu'il existe $\theta$ dans l'ensemble

$$
\Theta = \{ \theta = (\alpha_j, m_j, \sigma_j^2)_{1 \le j \le J} : \alpha_1 > 0, \dots, \alpha_J > 0, \alpha_1 + \dots + \alpha_J = 1 \}
$$

tel que les mesures $x_1, \dots, x_n$ de la taille de $n$ nids soient la réalisation d'un $n$-échantillon $X_1, \dots, X_n$ de loi de densité au point $x$ de la loi de $X$ par rapport à la mesure de Lebesgue $d\lambda$ de $\mathbb{R}$

$$
f_{\theta}(x) = \sum_{j=1}^{J} \alpha(j) \gamma_{m(j), v(j)}(x).
$$

On note $\bar{X}$ le vecteur $(X_1, \dots, X_n)$ et $\bar{Z}$ le vecteur $(Z_1, \dots, Z_n)$. On note $dN$
la mesure de comptage sur $\mathbb{N}$. On introduit les notations suivantes :

— densité au point $(x, z)$ de la loi du couple $(X, Z)$ par rapport à $d\lambda \otimes dN$ :

$$
h_{\theta}(x, z) = \alpha(z) \gamma_{m(z), v(z)}(x), \quad x \in \mathbb{R}, z \in \{1, \dots, J\},
$$

— densité au point $x$ de la loi de $X$ par rapport à $d\lambda$ :

$$
f_{\theta}(x) = \sum_{j=1}^{J} \alpha(j) \gamma_{m(j), v(j)}(x), \quad x \in \mathbb{R},
$$

— densité au point $z$ de la loi de $Z$ par rapport à la mesure $dN$ :

$$
g_{\theta}(z) = \alpha(z), \quad z \in \{1, \dots, J\},
$$

<footnote n="1">1. «Finite mixture» en anglais.</footnote>
— densité au point $x$ de la loi de $X$ sachant $Z$ par rapport à $d\lambda$ :

$$f_{\theta}(x|Z = j) = \gamma_{m(j),v(j)}(x), \quad x \in \mathbb{R},$$

— densité au point $z$ de la loi de $Z$ sachant $X$ par rapport à $dN$ :

$$g_{\theta}(z|X = x) = \frac{\alpha(z)\gamma_{m(z),v(z)}(x)}{\sum_{j=1}^{J} \alpha(j)\gamma_{m(j),v(j)}(x)}, \quad z \in \{1, \dots, J\}.$$

La dernière expression s'obtient grâce à la définition d'une probabilité conditionnelle : en effet, pour tout $A$ borélien et tout $1 \le z \le J$,

$$\mathbb{P}(X \in A, Z = z) = \mathbb{P}(Z = z)\mathbb{P}(X \in A|Z = z) = \alpha(z) \int_{A} \gamma_{m(z),v(z)}(x) dx$$

et

$$\begin{align*}
\mathbb{P}(X \in A, Z = z) &= \mathbb{P}(Z = z | X \in A)\mathbb{P}(X \in A) \\
&= \mathbb{P}(Z = z | X \in A) \sum_{j=1}^{J} \alpha(j) \int_{A} \gamma_{m(j),v(j)}(x) dx,
\end{align*}$$

de sorte que si $A$ est de la forme $[x - \varepsilon, x + \varepsilon]$ avec $\varepsilon$ tendant vers 0, alors on obtient bien

$$\alpha(z)\gamma_{m(z),v(z)}(x) = \mathbb{P}(Z = z | X = x) \sum_{j=1}^{J} \alpha(j)\gamma_{m(j),v(j)}(x).$$

## 8.2 Situation favorable mais irréaliste

On suppose dans cette section que l'on observe à la fois $X$ et $Z$, c'est-à-dire à la fois la taille du nid et l'espèce qui a fait le nid. Estimer les paramètres inconnus est alors aisé grâce à la méthode du maximum de vraisemblance. La log-vraisemblance du modèle complet (logarithme de la densité de la loi de l'échantillon $(\bar{X}, \bar{Z})$ par rapport à la mesure $(d\lambda \otimes dN)^{\otimes n})$ s'écrit

$$\begin{align*}
L(\theta, \bar{X}, \bar{Z}) &= \log \prod_{i=1}^{n} h_{\theta}(X_i, Z_i) \\
&= \sum_{i=1}^{n} [\log \alpha(Z_i) + \log \gamma_{m(Z_i),v(Z_i)}(X_i)].
\end{align*}$$

Pour un échantillon de taille $n$, notons, pour tout $1 \le j \le J$,

$$A_j = \{i = 1, \dots, n, Z_i = j\} \text{ et } C_j = \text{card}(A_j).$$

En regroupant les termes en fonction des valeurs possibles de $\bar{Z}$, on obtient

$$L(\theta, \bar{X}, \bar{Z}) = \sum_{j=1}^{J} C_j \log \alpha(j) + \sum_{j=1}^{J} \sum_{i \in A_j} \log \gamma_{m(j),v(j)}(X_i).$$
**Théorème 8.1** (Estimation pour modèle complet). La vraisemblance complète atteint son maximum en un point unique $\hat{\theta} \in \Theta$ donné par :

$$ \hat{\alpha}(j) = \frac{C_j}{n}, \quad \hat{m}(j) = \frac{1}{C_j} \sum_{i \in A_j} X_i \quad \text{et} \quad \hat{v}(j) = \frac{1}{C_j} \sum_{i \in A_j} (X_i - \hat{m}(j))^2. $$

*Démonstration.* Sous la contrainte $\alpha(1) + \dots + \alpha(J) = 1$, la fonction $\theta \mapsto L(\theta, \bar{X}, \bar{Z})$ admet un unique point critique, qui s'avère être un maximum. $\square$

Ce résultat ne répond pas à la question initiale puisque, en pratique, ne sont observées que les tailles des nids $x_1, \dots, x_n$. Il faudrait plutôt étudier la log-vraisemblance de l'échantillon $X_1, \dots, X_n$, notée $L_{\text{obs}}$ pour log-vraisemblance des observations, qui s'écrit :

$$ L_{\text{obs}}(\theta, \bar{X}) = \log \prod_{i=1}^{n} f_{\theta}(X_i) = \sum_{i=1}^{n} \log \left[ \sum_{j=1}^{J} \alpha(j) \gamma_{m(j),v(j)}(X_i) \right]. $$

Trouver un jeu de paramètres $\theta$ qui maximise cette quantité n'est pas facile.

## 8.3 Résolution du vrai problème

Ne disposant pas de la log-vraisemblance complète $L$, on la remplace par son espérance conditionnelle sachant les observations : on définit la log-vraisemblance conditionnelle des observations sous la loi de paramètre $\tilde{\theta}$, que nous noterons $L_c(\theta, \tilde{\theta}, \bar{X})$, par :

$$ L_c(\theta, \tilde{\theta}, \bar{X}) = \mathbb{E}_{\tilde{\theta}}(L(\theta, \bar{X}, \bar{Z})|\bar{X}) = \sum_{i=1}^{n} \sum_{j=1}^{J} \log h_{\theta}(X_i, j) g_{\tilde{\theta}}(j|X = X_i). $$

Voici l'expression de la vraisemblance conditionnelle $L_c$ dans notre cadre de mélange gaussien.

**Théorème 8.2** (Log-vraisemblance conditionnelle). La fonction $L_c$ est de la forme suivante :

$$ L_c(\theta, \tilde{\theta}, \bar{X}) = -\frac{n}{2} \log(2\pi) + \sum_{j=1}^{J} \left( \sum_{i=1}^{n} g_{\tilde{\theta}}(j|X = X_i) \right) \log \alpha(j) \\ - \frac{1}{2} \sum_{j=1}^{J} \left( \sum_{i=1}^{n} g_{\tilde{\theta}}(j|X = X_i) \right) \left( \log v(j) + \frac{(X_i - m(j))^2}{v(j)} \right). $$

Il est aisé de déterminer le maximum de la vraisemblance conditionnelle.
**Théorème 8.3** (Maximum de la log-vraisemblance conditionnelle). La fonction $\theta \mapsto L_c(\theta, \tilde{\theta}, \bar{X})$ admet un unique maximum $\theta_M$ donné par

$$ 
\begin{align*}
\alpha_M(j) &= \frac{1}{n} \sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i) \\
m_M(j) &= \frac{\sum_{i=1}^{n} X_i g_{\tilde{\theta}}(j | X = X_i)}{\sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i)} \\
v_M(j) &= \frac{\sum_{i=1}^{n} (X_i - m_M(j))^2 g_{\tilde{\theta}}(j | X = X_i)}{\sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i)}.
\end{align*}
$$

*Démonstration.* Sous la contrainte $\alpha(1) + \cdots + \alpha(J) = 1$, on trouve

$$ \alpha_M(j) = \frac{\sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i)}{\sum_{l=1}^{J} \sum_{i=1}^{n} g_{\tilde{\theta}}(l | X = X_i)} = \frac{1}{n} \sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i). $$

D'autre part $(m_M, v_M)$ est un point critique du dernier terme de la log-vraisemblance conditionnelle. $\square$

L'algorithme EM consiste à répéter successivement deux étapes consécutives comme suit :

- étape E(xpectation) : étant donnée une valeur $\theta_k$ du paramètre, on calcule la log-vraisemblance conditionnelle des observations $L_c(\theta, \theta_k, \bar{X})$ avec le théorème 8.2,

- étape M(aximization) : grâce au théorème 8.3, on choisit $\theta_{k+1}$ pour que la fonction $\theta \mapsto L_c(\theta, \theta_k, \bar{X})$ soit maximale au point $\theta_{k+1}$.

En pratique, étant données les observations $x_1, \dots, x_n$ et des valeurs initiales des paramètres, l'algorithme consiste à répéter les calculs suivants :

- à partir du paramètre $\theta_k$, on calcule la matrice $H^{(k)}$ de taille $n \times J$ suivante

$$ H_{ij}^{(k)} = \frac{\alpha_k(j) \gamma_{m_k(j), v_k(j)}(X_i)}{\sum_{l=1}^{J} \alpha_k(l) \gamma_{m_k(l), v_k(l)}(X_i)}, $$

- on en déduit $\theta_{k+1}$ par les relations suivantes :

$$ 
\begin{align*}
\alpha_{k+1}(j) &= \frac{1}{n} \sum_{i=1}^{n} H_{ij}^{(k)}, \\
m_{k+1}(j) &= \frac{\sum_{i=1}^{n} X_i H_{ij}^{(k)}}{\sum_{i=1}^{n} H_{ij}^{(k)}}, \\
v_{k+1}(j) &= \frac{\sum_{i=1}^{n} (X_i - m_{k+1}(j))^2 H_{ij}^{(k)}}{\sum_{i=1}^{n} H_{ij}^{(k)}}.
\end{align*}
$$

Le théorème 8.3 ne vaut que si nous savons estimer $\tilde{\theta}$. Dans l'algorithme EM, ceci est fait numériquement. Le résultat suivant montre que la log-vraisemblance $L_{\text{obs}}$ est croissante le long de l'algorithme.
**Théorème 8.4** (Croissance de la vraisemblance). La suite $(\theta_k)_{k \ge 0}$ construite par l'algorithme EM vérifie la propriété de stabilité numérique suivante :

$$L_{\text{obs}}(\theta_{k+1}, \bar{X}) \ge L_{\text{obs}}(\theta_k, \bar{X}),$$

où $L_{\text{obs}}(\theta, \bar{X})$ est la log-vraisemblance des observations.

**Démonstration.** Posons

$$H(\theta, \theta_k) = \mathbb{E}_{\theta_k}(\log g_{\theta}(\bar{Z} | \bar{X})) = \sum_{i=1}^{n} \mathbb{E}_{\theta_k}(\log g_{\theta}(Z_i | X_i)).$$

Puisque $h_{\theta}(x, z) = f_{\theta}(x)g_{\theta}(z|X=x)$, la vraisemblance conditionnelle s'écrit

$$L_c(\theta, \theta_k, \bar{X}) = L_{\text{obs}}(\theta, \bar{X}) + H(\theta, \theta_k).$$

Lors de l'étape M, la fonction $\theta \mapsto L_c(\theta, \theta_k, \bar{X})$ est maximisée en $\theta$, d'où

$$L_c(\theta_{k+1}, \theta_k, \bar{X}) \ge L_c(\theta_k, \theta_k, \bar{X}).$$

On en déduit que

$$L_{\text{obs}}(\theta_{k+1}, \bar{X}) - L_{\text{obs}}(\theta_k, \bar{X}) \ge H(\theta_k, \theta_k) - H(\theta_{k+1}, \theta_k).$$

Par définition de $H$,

$$H(\theta_k, \theta_k) - H(\theta_{k+1}, \theta_k) = \mathbb{E}_{\theta_k} \log \frac{g_{\theta_k}(\bar{Z} | \bar{X})}{g_{\theta_{k+1}}(\bar{Z} | \bar{X})}.$$

En notant $\mu_{\theta}$ la loi de $\bar{Z}$ sachant $\bar{X}$ sous $\mathbb{P}_{\theta}$, on remarque que la quantité ci-dessus n'est rien d'autre que l'entropie relative de $\mu_{\theta_k}$ par rapport à $\mu_{\theta_{k+1}}$ :

$$H(\theta_k, \theta_k) - H(\theta_{k+1}, \theta_k) = \int \log \frac{d\mu_{\theta_k}}{d\mu_{\theta_{k+1}}} d\mu_{\theta_k} = \int \log \frac{d\mu_{\theta_k}}{d\mu_{\theta_{k+1}}} \log \frac{d\mu_{\theta_k}}{d\mu_{\theta_{k+1}}} d\mu_{\theta_{k+1}}.$$

L'inégalité de Jensen assure que cette quantité est positive. $\square$

Le théorème 8.4 montre que la suite $(L_{\text{obs}}(\theta_k, \bar{X}))_{k \ge 0}$ est croissante. Elle converge donc vers un point critique (maximum local ou point selle) de la vraisemblance $L_{\text{obs}}(\cdot, \bar{X})$. Contrairement au cas où l'on observe en même temps $X$ et $Z$, il peut exister des maxima locaux qui vont piéger l'algorithme, ce qui rend le résultat sensible au point de départ. On peut améliorer le comportement de l'algorithme en incorporant un peu d'aléa, comme par exemple dans l'algorithme du recuit simulé du chapitre 4, mais cela est un peu illuisoire. De plus on ne sait pas combien d'itérations sont requises pour approcher suffisamment un maximum local. Il n'y a pas de critère d'arrêt performant.

La figure 8.1 donne l'estimation du mélange de trois gaussiennes par l'algorithme EM avec 1000 itérations et les paramètres

$$\alpha = (0.4 \ 0.4 \ 0.2), \ m = (-2 \ 0 \ 2), \ v = (0.3 \ 0.2 \ 0.2), \ n = 1000.$$
<figure><img src="image_26.png" /><figcaption>Fig. 8.1. Estimation du mélange de trois gaussiennes par l'algorithme EM.</figcaption></figure>

## 8.4 Pour aller plus loin

Le théorème 8.3 se généralise aux familles exponentielles, de vraisemblance totale

$$L(\theta, \bar{x}, \bar{z}) = \sum_{i=1}^{n} \log h_{\theta}(x_i, z_i) = \sum_{i=1}^{n} \log g_{\theta}(z_i) + \sum_{i=1}^{n} \log f_{\theta}(x_i | Z = z_i).$$

Après intégration, il vient en rappelant que $g_{\theta}(z) = \alpha(z)$,

$$L_c(\theta, \tilde{\theta}, \bar{x}) = \sum_{j=1}^{J} \left( \sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = x_i) \right) \log \alpha(j) + \sum_{j=1}^{J} \sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = x_i) \log f_{\theta}(x_i | Z = j).$$

Si $\theta = (\alpha, \beta)$ où $\alpha$ est encore la loi mélange et $\beta(j)$ représente les paramètres de la loi $\text{Loi}(X | Z = j)$, on détermine séparément $\alpha_M$ (qui a toujours la même expression) et $\beta_M$ qui s'obtient quasiment comme pour le maximum de la vraisemblance totale : on a seulement pondéré la contribution de l'observation $X_i$ par le coefficient $g_{\tilde{\theta}}(j | X = X_i)$. Par exemple, si les composantes du mélange $(\text{Loi}(X | Z = j))_{1 \le j \le J}$ sont des lois exponentielles de paramètres $(\lambda(i))_{1 \le i \le J}$, la log-vraisemblance conditionnelle s'écrit
$$
L_c(\theta, \tilde{\theta}, \bar{X}) = \sum_{j=1}^{J} \left( \sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i) \right) \log \alpha(j) \\
+ \sum_{j=1}^{J} \sum_{i=1}^{n} g_{\tilde{\theta}}(j | X = X_i) (\log \lambda(j) - \lambda(j) x_i).
$$

Un pas de l’algorithme EM a la forme suivante :

$$
H_{ij}^{(k)} = \frac{\alpha_k(j) \lambda_k(j) e^{-\lambda_k(j)x_i}}{\sum_{l=1}^{J} \alpha_k(l) \lambda_k(l) e^{-\lambda_k(l)x_i}},
$$
$$
\alpha_{k+1}(j) = \frac{1}{n} \sum_{i=1}^{n} H_{ij}^{(k)},
$$
$$
\lambda_{k+1}(j) = \frac{\sum_{i=1}^{n} H_{ij}^{(k)}}{\sum_{i=1}^{n} X_i H_{ij}^{(k)}}.
$$

Si de nombreux travaux existaient déjà autour de cette question, c'est vérita-
blement l'article [DLR77] de Arthur Dempster, Nan Laird, et Donald Rubin
qui définit pour la première fois l'algorithme EM dans un cadre général. On
trouvera dans la bibliographie de ce travail les références aux résultats an-
térieurs. Depuis, cet algorithme est très souvent utilisé dans des cadres et
sous des formes variés. Pour réduire la dépendance aux conditions initiales,
de nombreuses versions randomisées de l'algorithme EM ont été développées,
inspirées notamment de l'algorithme d'approximation stochastique de Kiefer-
Wolfowitz pour le calcul du maximum de fonction sous forme d'espérance,
lui même inspiré de l'algorithme d'approximation stochastique de Robbins-
Monro pour le calcul de zéro de fonction sous forme d'espérance. On pourra
par exemple consulter à ce sujet l'article [DLM99] de Bernard Delyon, Marc
Lavielle, et Éric Moulines, et le livre [Duf97] de Marie Duflo. Lorsqu'il n'est pas
possible de calculer la log-vraisemblance conditionnelle, on peut utiliser une
méthode de Monte-Carlo comme l'algorithme de Metropolis-Hastings du cha-
pitre 5 pour obtenir une approximation de cette fonction dont on cherche en-
suite un maximum numériquement. Ce type d'approche est étudié par exemple
par Estelle Kuhn et Marc Lavielle dans [KL04].

Une étude de la dégénérescence de la vraisemblance des mélanges finis
gaussiens et ses liens avec l'algorithme EM est menée par Christophe Biernacki
et Stéphane Chrétien dans [BC03].
# Urnes d'Ehrenfest

**Mots-clés.** Urne; mesure invariante; mesure réversible; temps de retour; convergence à l'équilibre; convergence abrupte.

**Outils.** Chaîne de Markov; chaîne de Markov à temps continu; couplage; distance en variation totale; décomposition spectrale; collectionneur de coupons; inégalité de Chernoff; inégalité de Tchebychev; théorème limite central.

## Difficulté. **

Le second principe de la thermodynamique postule que pour tout système isolé, une grandeur macroscopique extensive appelée *entropie* augmente au cours du temps. Dans la théorie cinétique des gaz de Boltzmann, le second principe est une conséquence d'une approche statistique de la mécanique classique appliquée aux constituants microscopiques. Bien que séduisante, cette approche fut critiquée par les adversaires de l'atomisme. En effet, comment lever la contradiction entre la *réversibilité microscopique* héritée de la mécanique classique et l'*irréversibilité macroscopique* de la thermodynamique ? Le modèle d'Ehrenfest propose, dans un cas extrêmement simple, une explication à ce paradoxe. L'outil clé ici est celui des chaînes de Markov finies.

## 9.1 Modèle microscopique

Considérons deux récipients *A* et *B*, connectés par un petit trou, et contenant un gaz constitué de *a* molécules (figure 9.1). À chaque instant, l'une des *a* molécules, au hasard, passe d'un récipient à l'autre en empruntant le petit trou. On souhaite décrire l'évolution du système. On modélise ces deux récipients par deux urnes *A* et *B* dans lesquelles sont réparties *a* boules numérotées de 1 à *a*. On associe à une configuration, c'est-à-dire à une répartition des *a* boules, un élément *x* = (*x*⟨1⟩, ..., *x*⟨*a*⟩) du cube discret
<figure><img src="image_27.png" /><figcaption>Fig. 9.1. Répartition des molécules dans les deux urnes.</figcaption></figure>

$$F := \{0, 1\}^a$$

où $x(i) = 1$ si la boule $i$ est dans l'urne $A$ et $x(i) = 0$ sinon. On dit que deux configurations $x$ et $y$ de $F$ sont voisines, et on note $x \sim y$ (ou encore $x \sim_i y$), si elles ne diffèrent que d'une coordonnée (qui est $i$). Cette relation de voisinage définit sur $F$ une structure de graphe non orienté dont les sommets sont les points de $F$ et les arêtes sont les paires de sommets voisins. La marche aléatoire sur le cube discret $Y = (Y_n)_{n \ge 0}$ est la chaîne de Markov d'espace d'états $F$, dont la matrice de transition est donnée, pour tous $x, y \in F$, par

$$Q(x, y) = \begin{cases} \frac{1}{a} & \text{si } x \sim y, \\ 0 & \text{sinon.} \end{cases}$$

Cette matrice est symétrique et en particulier doublement stochastique. Le temps d'atteinte de $y \in E$ est la variable aléatoire $\inf\{n \ge 0 : Y_n = y\}$. Sur l'événement $\{Y_0 = y\}$ le temps de retour en $y$ est $\inf\{n \ge 1 : Y_n = y\}$.

**Lemme 9.1** (Premières propriétés). La chaîne de Markov Y d'espace d'états $F$ et de matrice de transition $Q$ est irréductible, récurrente, de période 2, et la loi uniforme sur $F$ notée $\pi$ est invariante et réversible<fnref n="1" />. De plus, le temps moyen de retour en $x \in F$ ne dépend pas de $x$ et vaut $2^a : \text{pour tout } x \in F$,

<footnote n="1">1. C'est-à-dire que $\pi(x)Q(x, y) = \pi(y)Q(y, x)$ pour tous $x, y$. Sur les trajectoires de $Y$, si $Y_0 \sim \pi$ alors $(Y_0, \dots, Y_n)$ et $(Y_n, \dots, Y_0)$ ont même loi, pour tout $n \in \mathbb{N}$.</footnote>
$$
\mathbb{E}(\inf\{n \ge 1 : Y_n = x\} | Y_0 = x) = \frac{1}{\pi(x)} = 2^a.
$$

La loi du temps d'atteinte de y en partant de x ne dépend que de

$$
d(x, y) = \sum_{i=1}^{a} |x_i - y_i| = \|x - y\|_{1}.
$$

Cette quantité est une distance qui compte le nombre de coordonnées diffé-
rentes entre x et y. On note $m_d$ (en oubliant un moment la dépendance en $a$)
le temps d'atteinte moyen de $y$ partant de $x$ lorsque $d(x, y) = d$, c'est-à-dire

$$
m_d := \mathbb{E}(\inf\{n \ge 0 : Y_n = y\} | Y_0 = x).
$$

**Théorème 9.2** (Temps moyens). *Pour* $1 \le d \le a$, on a

$$
m_d = \sum_{i=1}^{d} \mathbf{Q}_{a-i}^{a} \quad \text{où} \quad \mathbf{Q}_{i}^{a} = \sum_{k=0}^{i} \frac{\binom{a}{k}}{\binom{a-1}{i}} \quad \text{si } 1 \le i \le a-1.
$$

Démonstration. D'après la propriété de Markov, la suite $(m_d)_{0 \le d \le a+1}$ vérifie

$$
\begin{cases}
m_d = 1 + \frac{d}{a}m_{d-1} + \frac{a-d}{a}m_{d+1} & \text{si } 1 \le d \le a, \\
m_0 = m_{a+1} = 0, &
\end{cases}
$$

et ce système d'équations de récurrence possède une unique solution. $\square$

Comme $Y$ est de période 2, la loi de $Y_n$ ne converge pas vers la loi invariante $\pi$ quand $n \to \infty$. Pour obtenir une chaîne apériodique, on la rend paresseuse $^2$. Cela consiste à modifier la matrice de transition en posant, pour tous $x, y \in F$

$$
\mathbf{R}(x, y) = \frac{1}{2} \mathbf{1}_{\{y=x\}} + \frac{1}{2} \mathbf{Q}(x, y).
$$

La chaîne $(Z_n)_{n \ge 0}$ associée à la matrice de transition $\mathbf{R}$ reste sur place avec
probabilité $1/2$ ou se déplace avec la même probabilité vers l'un des sommets
voisins choisi selon la loi uniforme. Cette chaîne est irréductible récurrente
apériodique et réversible pour la loi $\pi$. La loi de $Z_n$ converge donc vers $\pi$
quelle que soit la loi de $Z_0$. Nous présentons deux méthodes pour quantifier
cette convergence. La première fait appel à une méthode de couplage et au
collectionneur de coupons de la section 1.5. La seconde, plus fine, utilise la
décomposition spectrale de $\mathbf{R}$. On note $d_{\text{VT}}(\cdot, \cdot)$ la distance en variation totale
étudiée dans la section 1.3.

**Théorème 9.3** (Temps long et collectionneur de coupons). *Pour tout* $n \ge 0$,

$$
\max_{x \in F} d_{\text{VT}} (\mathbf{R}^n(x, \cdot), \pi) \le P(T_a > n),
$$

où $T_a$ est le temps de complétion d'une collection de a coupons.

<footnote n="2">2. « *Lazy chain* » en anglais.</footnote>
Démonstration. Fixons $n \ge 0$. Grâce au théorème 1.5, pour tout $x \in F$,

$$
\begin{align*}
d_{\mathrm{VT}}(\mathbf{R}^n(x, \cdot), \pi) &= d_{\mathrm{VT}}(\mathbf{R}^n(x, \cdot), \pi \mathbf{R}^n) \\
&= \frac{1}{2} \sum_{y \in F} \left| \mathbf{R}^n(x, y) - \sum_{x' \in F} \pi(x') \mathbf{R}^n(x', y) \right| \\
&\le \frac{1}{2} \sum_{y \in F} \sum_{x' \in F} \pi(x') \left| \mathbf{R}^n(x, y) - \mathbf{R}^n(x', y) \right| \\
&\le \max_{x' \in F} d_{\mathrm{VT}}(\mathbf{R}^n(x, \cdot), \mathbf{R}^n(x', \cdot)).
\end{align*}
$$

Pour contrôler $d_{\mathrm{VT}}(\mathbf{R}^n(x, \cdot), \mathbf{R}^n(x', \cdot))$, on construit un couple $(Z, Z')$ de chaînes de Markov $Z := (Z_n)_{n \ge 0}$ et $Z' := (Z'_n)_{n \ge 0}$ de même matrice de transition $\mathbf{R}$ et de conditions initiales $Z_0 = x$ et $Z'_0 = x'$, et dont les trajectoires sont égales après un temps aléatoire (de coalescence) que l'on sait contrôler. Plus précisément on se donne deux suites indépendantes de variables aléatoires indépendantes $(U_n)_{n \ge 0}$ et $(V_n)_{n \ge 0}$ de lois respectives de Bernoulli $\mathrm{Ber}(1/2)$ et uniforme sur $\{1, \dots, a\}$. On pose alors pour tout $n \ge 0$ et $1 \le i \le a$,

$$
\begin{cases}
Z_{n+1}(V_n) = U_n, \\
Z_{n+1}(i) = Z_n(i) \text{ si } i \neq V_n,
\end{cases}
\quad
\begin{cases}
Z'_{n+1}(V_n) = U_n, \\
Z'_{n+1}(i) = Z'_n(i) \text{ si } i \neq V_n.
\end{cases}
$$

On vérifie aisément que $Z$ et $Z'$ sont bien des chaînes de Markov de matrice
de transition $\mathbf{R}$. À chaque instant une coordonnée est choisie uniformément
et elle est positionnée à une même valeur aléatoire pour les deux chaînes. Si
les coordonnées étaient égales, elles le restent. Sinon, elles le deviennent. La
loi du temps de coalescence global

$$
T_d := \inf\{n \ge 0 : Z_n = Z'_n\}
$$

ne dépend que du nombre $d := d(x, x')$ de coordonnées différentes entre $x$ et $x'$. Pour toute fonction $f : F \to \mathbb{R}$, il vient, en notant $\Delta_n := f(Z_n) - f(Z'_n)$,

$$
\mathrm{E}(\Delta_n) = \underbrace{\mathrm{E}(\Delta_n \mathbf{1}_{\{T_d \le n\}})}_{=0} + \mathrm{E}(\Delta_n \mathbf{1}_{\{T_d > n\}}) = \mathrm{E}(\Delta_n \mathbf{1}_{\{T_d > n\}}),
$$

ce qui donne $\mathrm{E}(\Delta_n) \le 2\|f\|_\infty \mathbb{P}(T_d > n)$, d'où, grâce au théorème 1.5,

$$
\begin{align*}
d_{\mathrm{VT}}(\mathbf{R}^n(x, \cdot), \mathbf{R}^n(x', \cdot)) &= d_{\mathrm{VT}}(\mathrm{Loi}(Z_n), \mathrm{Loi}(Z'_n)) \\
&= \frac{1}{2} \sup_{\|f\|_{\infty} \le 1} \mathrm{E}(f(Z_n) - f(Z'_n)) \\
&\le \mathbb{P}(T_d > n).
\end{align*}
$$

Alternativement, on peut utiliser le théorème 1.9 et la définition de $T_d$ :

$$
d_{\mathrm{VT}}(\mathrm{Loi}(Z_n), \mathrm{Loi}(Z'_n)) \le \mathbb{P}(Z_n \ne Z'_n) \le \mathbb{P}(T_d > n).
$$
À présent, pour contrôler la queue de distribution de $T_d$, on observe tout d'abord que $T_d$ coïncide avec le premier instant où toutes les coordonnées dans $I := \{i \in \{1, \dots, a\} : x_i \neq x'_i\}$ ont été mises à jour au moins une fois. On a card($I$) = $d$. Supposons pour simplifier la notation que $I = \{1, \dots, d\}$. Si les variables aléatoires $(\tau_i)_{1 \le i \le d}$ désignent les temps successifs de modification d'une nouvelle coordonnée alors $T_d = \tau_d$ et les variables aléatoires $(G_i)_{1 \le i \le d}$ définies par $G_1 = 1$ et $G_i = \tau_i - \tau_{i-1}$ pour $2 \le i \le d$ sont indépendantes et de lois géométriques de paramètres respectifs $(d+1-i)/d$. On retrouve ici le collectionneur de coupons de la section 1.5. Il ne reste plus qu'à observer que $\mathbb{P}(T_d > n)$ croît en $d$, et atteint son maximum pour $d = a$, ce qui revient à prendre $x = (0, \dots, 0)$ et $x' = (1, \dots, 1)$. $\square$

Le théorème 1.14 assure que pour tout réel $c \ge 0$,

$$ \mathbb{P}(T_a \ge a \log(a) + ca) \le e^{-c+1/a} $$

Le théorème suivant raffine cette estimation de $\max_{x \in F} d_{VT}(\mathbf{R}^n(x, \cdot), \pi)$.

**Théorème 9.4 (Temps long).** Si

$$ n = \frac{a \log(a)}{2} + ca \quad \text{avec } c \ge 0 $$

alors

$$ \max_{x \in F} d_{VT}(\mathbf{R}^n(x, \cdot), \pi) \le \frac{e^{-c}}{\sqrt{2}} $$

*Démonstration.* On munit $\mathbb{R}^F \equiv \mathbb{R}^{2^a}$ du produit scalaire usuel noté $\langle \cdot, \cdot \rangle$. La matrice stochastique $\mathbf{R}$ est symétrique$^3$, son spectre est réel et inclus dans

$$ [-1, 1] = \mathbb{R} \cap \{z \in \mathbb{C} : |z| \le 1\}. $$

De plus, il existe une base orthonormée de $\mathbb{R}^F$ formée de vecteurs propres $(f_j)_{1 \le j \le 2^a}$ de $\mathbf{R}$ associés aux valeurs propres $(\lambda_j)_{1 \le j \le 2^a}$ classées par ordre décroissant. En particulier, $\lambda_1$ vaut 1 et $f_1$ est la fonction constante égale à $2^{-a/2}$. Rappelons que $\pi$ est la loi uniforme sur $F$ de cardinal $2^a$ et donc $\pi(x) = 2^{-a}$ pour tout $x \in F$. Pour toute fonction $f$ de $F$ dans $\mathbb{R}$ et $x$ de $F$,

$$ \mathbf{R}^n(f)(x) = \sum_{j=1}^{2^a} \langle f, f_j \rangle f_j(x) \lambda_j^n = \pi(f) + \sum_{j=2}^{2^a} \langle f, f_j \rangle f_j(x) \lambda_j^n. $$

Pour la fonction $z \in F \mapsto f(z) = \mathbb{1}_y(z)$, on obtient la décomposition spectrale

<footnote n="3">3. Pour toute matrice de transition sur un espace d'états fini, la loi uniforme est invariante (respectivement réversible) si et seulement si la matrice est doublement stochastique (respectivement symétrique).</footnote>
$$ \mathbf{R}^n(x, y) = \frac{1}{2^a} + \sum_{j=2}^{2^a} f_j(x) f_j(y) \lambda_j^n. $$

À présent, d'après le théorème 1.5 et l'inégalité de Schwarz,

$$ 4d_{\text{VT}}(\mathbf{R}^n(x, \cdot), \pi(\cdot))^2 = \left(\sum_{y \in F} |\mathbf{R}^n(x, y) - \pi(y)|\right)^2 \le 2^a \sum_{y \in F} |\mathbf{R}^n(x, y) - \pi(y)|^2. $$

Puisque $(f_i)_i$ est une base orthonormée, la décomposition spectrale donne

$$ \sum_{y \in F} |\mathbf{R}^n(x, y) - \pi(y)|^2 = \sum_{i,j=2}^{2^a} f_i(x) f_j(x) \langle f_i, f_j \rangle \lambda_i^n \lambda_j^n = \sum_{i=2}^{2^a} f_i(x)^2 \lambda_i^{2n}. $$

Par symétrie de $F$ et $\mathbf{R}$, la quantité $\sum_{y \in F} |\mathbf{R}^n(x, y) - \pi(y)|^2$ ne dépend pas de $x \in F$. En sommant sur $x \in F$, on obtient

$$ 2^a \sum_{y \in F} |\mathbf{R}^n(x, y) - \pi(y)|^2 = \sum_{i=2}^{2^a} \lambda_i^{2n}, $$

qui fournit la majoration

$$ 4d_{\text{VT}}(\mathbf{R}^n(x, \cdot), \pi(\cdot))^2 \le \sum_{i=2}^{2^a} \lambda_i^{2n}. $$

Or on peut établir que les valeurs propres de $\mathbf{R}$ sont les réels $(1 - k/a)_{0 \le k \le a}$ avec les multiplicités respectives $\binom{a}{k}$ pour $0 \le k \le a$. On a donc

$$ \begin{aligned} 4d_{\text{VT}}(\mathbf{R}^n(x, \cdot), \pi)^2 &\le \sum_{k=1}^a \left(1 - \frac{k}{a}\right)^{2n} \binom{a}{k} \\ &\le \sum_{k=1}^a e^{-2nk/a} \binom{a}{k} = \left(1 + e^{-2n/a}\right)^a - 1. \end{aligned} $$

Le choix $n = (1/2)a \log(a) + ca$ fournit

$$ 4d_{\text{VT}}(\mathbf{R}^n(x, \cdot), \pi)^2 \le \left(1 + \frac{e^{-2c}}{a}\right)^a - 1 \le e^{e^{-2c}} - 1. $$

On a $e^{e^{-2c}} - 1 \le 2e^{-2c}$ dès que $c \ge 0$, d'où le résultat. $\square$

<figure><img src="image_28.png" /></figure>

Une autre façon de contourner le problème de la périodicité de la matrice $\mathbf{Q}$ est d'introduire une version à temps continu $(Z_t)_{t \ge 0}$ de la chaîne $(Y_n)_{n \ge 0}$. On suppose que chaque coordonnée subit une inversion 4 de 0 à 1 ou de 1 à 0 à

<footnote n="4">4. «Flip» en anglais.</footnote>
des instants aléatoires. Plus précisément, soient $(N^{(i)})_{1\le i\le a}$ des processus de Poisson indépendants de même intensité $1/a$. Posons

$$Z_t(i) = Z_0(i) + N_t^{(i)} \mod 2.$$

Chaque coordonnée est modifiée indépendamment des autres à des temps exponentiels de paramètre $1/a$. Cette dynamique peut être résumée par le générateur infinitésimal du processus $\mathbf{L} = \mathbf{Q} - \mathbf{I}$ : pour toute fonction $f$ de $F$ dans $\mathbb{R}$ et tout $x \in F$,

$$\mathbf{L}f(x) = \frac{1}{a} \sum_{y \sim x} (f(y) - f(x)) = \frac{1}{a} \sum_{i=1}^{n} (f(x_{(i)}) - f(x))$$

où $x_{(i)} \in F$ s'obtient à partir de $x$ en inversant la $i$-e coordonnée. La loi uniforme $\pi$ est solution de $\pi \mathbf{L} = 0$ et constitue la loi invariante de $Z$.

**Théorème 9.5** (Temps long du modèle à temps continu). Pour tout $t \ge 0$,

$$\max_{x \in F} d_{\text{VT}} (\text{Loi}(Z_t | Z_0 = x), \pi) \le \mathbb{P}(S_a > t)$$

où $S_a := \max_{1 \le i \le a} E_i$ et $E_1, \dots, E_a$ sont des v.a.r. i.i.d. de loi Exp$(2/a)$.

**Démonstration.** On reprend la méthode de preuve du théorème 9.3 par couplage. On montre tout d'abord de la même manière que pour tout $x \in F$,

$$d_{\text{VT}}(\text{Loi}(Z_t | Z_0 = x), \pi) \le \max_{x' \in F} d_{\text{VT}}(\text{Loi}(Z_t | Z_0 = x), \text{Loi}(Z_t | Z_0 = x')).$$

Ensuite on construit un couplage comme suit : les couples de coordonnées $(Z.(i), Z'(i))_{1\le i\le a}$ évoluent indépendamment et pour chaque $i$, $Z_t(i)$ et $Z'_t(i)$ sont indépendants tant qu'ils sont différents et restent collés lorsqu'ils se touchent. Le générateur infinitésimal $\mathbb{L}$ du processus $((Z_t, Z'_t))_{t\ge 0}$ est donné pour toute fonction $f$ de $F \times F$ dans $\mathbb{R}$ et tout $(x, x') \in F \times F$ par

$$\mathbb{L}f(x, x') = \frac{1}{a} \sum_{i=1}^{a} \left( \begin{aligned} & (f(x_{(i)}, x'_{(i)}) - f(x, x')) \mathbb{1}_{x_i=x_i'} \\ & + (f(x, x'_{(i)}) - f(x, x')) \mathbb{1}_{x_i \neq x_i'} \\ & + (f(x_{(i)}, x') - f(x, x')) \mathbb{1}_{x_i \neq x_i'} \end{aligned} \right).$$

Il coïncide avec $\mathbf{L}$ lorsque $f$ ne dépend que de la première ou que de la seconde variable : $Z$ et $Z'$ sont des chaînes de générateur $\mathbf{L}$. On prend pour condition initiale $(Z_0, Z'_0) = (x, x')$. Comme dans la preuve du théorème 9.3 par couplage, la loi du temps de coalescence global

$$S_d := \inf\{t \ge 0 : Z_t = Z'_t\}$$

ne dépend de $(x, x')$ qu'à travers $d := d(x, x')$, et
$$d_{VT}(\text{Loi}(Z_t | Z_0 = x), \text{Loi}(Z_t | Z_0 = x')) \le P(S_d > t).$$

Pour contrôler la queue de distribution de $S_d$, on recherche une représentation plus commode. Il y a au départ $a-d+2d=a+d$ transitions possibles, de même taux $1/a$, dont $2d$ provoquent une nouvelle coalescence de coordonnée. Cela donne un taux total de $\lambda = (a+d)/a$ et une probabilité de coalescence de coordonnée de $p = 2d/(a+d)$. La construction des trajectoires<fnref n="5" /> et une probabilité des lois exponentielles<fnref n="6" /> font que la première coalescence de coordonnée a lieu au bout d'une durée aléatoire $F_d$ de loi exponentielle de paramètre

$$p\lambda = \frac{2d}{a+d} \frac{a+d}{a} = \frac{2}{a}d.$$

Plus généralement, pour tout $i \in \{1, \dots, d\}$, la $i^e$ coalescence de coordonnée a lieu au bout d'une durée aléatoire $F_{d-i+1}$ de loi exponentielle de paramètre

$$\frac{2(d-i+1)}{a+d-i+1} \frac{a+d-i+1}{a} = \frac{2}{a}(d-i+1).$$

Ces durées sont indépendantes, et le temps de coalescence global

$$S_d := \inf\{t \ge 0 : Z_t = Z'_t\}$$

a la même loi que $F_d+\dots+F_1$ où $F_d, \dots, F_1$ sont des v.a.r. indépendantes avec $F_k \sim \text{Exp}(2k/a)$ pour tout $k \in \{1, \dots, d\}$. Le lemme de Rényi 11.3 indique que $F_d+\dots+F_1$ a la même loi que $\max(E_1, \dots, E_d)$ où $E_1, \dots, E_d$ sont des v.a.r. i.i.d. de loi exponentielle de paramètre $2/a$.

$$S_d \stackrel{\text{loi}}{=} \max_{1 \le i \le d} E_i.$$

Pour tout $t \ge 0$, la quantité $P(S_d > t)$ croît en $d$, et atteint son maximum $P(S_a > t)$ pour $d=a$, ce qui correspond à $x = (0, \dots, 0)$ et $x' = (1, \dots, 1)$. $\square$

Le temps aléatoire $S_a$ est plus facile à étudier que le temps $T_a$ de complétion d'une collection même s'ils ont le même comportement asymptotique comme le montre le résultat suivant.

**Théorème 9.6** (Temps long pour la chaîne à temps continu).

$$\frac{S_a - (a/2) \log(a)}{a/2} \xrightarrow{a \to \infty} \text{Gumbel}$$

où la loi de Gumbel a pour fonction de répartition $t \in \mathbb{R} \mapsto e^{-e^{-t}}$.

<footnote n="5">5. Si $Z$ est une chaîne de Markov à temps continu et à espace d'états discret $E$ de générateur infinitésimal $G = (G(x, y))_{x,y \in E}$ alors sachant que la chaîne est en $x$ au temps $t$, le prochain saut a lieu au bout d'un temps exponentiel de paramètre $-G(x, x)$ et se fait en $y \neq x$ avec probabilité $-G(x, y)/G(x, x)$.</footnote>

<footnote n="6">6. Le lemme 11.4 affirme que si $G, E_1, E_2, \dots$ sont des v.a. indépendantes avec $G$ de loi géométrique de paramètre $p$ et $E_1, E_2, \dots$ de loi exponentielle de paramètre $\lambda$ alors la somme aléatoire $E_1 + \dots + E_G$ suit la loi exponentielle de paramètre $p\lambda$.</footnote>
Tout comme pour $T_a$ (théorème 1.17) il y a convergence abrupte de $S_a$ autour de la valeur $(a/2) \log(a)$ avec une fenêtre de largeur $a$.

*Démonstration.* Pour tout $a \ge 1$ et $r \ge -\log(a)$,

$$ 
\begin{align*}
\mathbb{P}\left(\frac{S_a - (a/2) \log(a)}{a/2} \le r\right) &= \mathbb{P}\left(\max_{1 \le i \le a}(E_i) \le \frac{ar}{2} + \frac{a \log(a)}{2}\right) \\
&= \mathbb{P}\left(E_1 \le \frac{a}{2}(r + \log(a))\right)^a \\
&= \left(1 - \frac{e^{-r}}{a}\right)^a.
\end{align*}
$$

Ainsi, pour tout $r \in \mathbb{R}$,

$$ \mathbb{P}\left(\frac{S_a - (a/2) \log a}{a/2} \le r\right) \underset{a \to \infty}{\longrightarrow} e^{-e^{-r}}, $$

ce qui est la convergence annoncée.

□

## 9.2 Urne d'Ehrenfest

Observer l'évolution de la chaîne Y demanderait de pouvoir déterminer la position de chaque molécule ce qui est physiquement impossible car $a$ est de l'ordre du nombre d'Avogadro<fnref n="7" />. La grandeur macroscopique qui nous intéresse, et que l'on peut mesurer, est la proportion ou le nombre de molécules situées dans l'urne A. Pour simplifier certaines expressions dans la suite, on sera amené à se placer parfois dans le cas où $a$ est pair. La lettre $b$ désignera alors toujours dans la suite l'entier $a/2$. À la chaîne de Markov Y qui décrit le système microscopiquement, on associe le processus $(X_n)_{n \ge 0}$ à valeurs dans

$$ E := \{0, 1, \dots, a\} $$

où

$$ X_n := \sum_{i=1}^{a} Y_n(i) $$

est la somme des coordonnées de $Y_n$, qui décrit le système macroscopiquement.

**Théorème 9.7 (Chaîne d'Ehrenfest).** Le processus $(X_n)_{n \ge 0}$ est une chaîne de Markov sur E de matrice de transition

<footnote n="7">7. Ce nombre approximativement égal à $6.022 \times 10^{23}$ définit la mole et correspond au nombre d'atomes dans 12 grammes de carbone 12. C'est aussi environ le nombre de molécules contenues dans 22.414 litres d'un gaz parfait dans les conditions normales de température (0° C) et de pression (1 Atmosphère).</footnote>
$$
\mathbf{P}(x, x') = \begin{cases} x/a & \text{si } x' = x - 1; \\ 1 - x/a & \text{si } x' = x + 1; \quad \text{pour tous } x, x' \in E, \\ 0 & \text{sinon,} \end{cases}
$$

irréductible, récurrente, de période 2. La loi binomiale Bin(a, 1/2) est inva-
riante et réversible.

Démonstration. Soit $\mathcal{F}_n$ et $\mathcal{G}_n$ les tribus engendrées par les suites $(X_k)_{0 \le k \le n}$
et $(Y_k)_{0 \le k \le n}$ respectivement, et soit $S(y) = y(1) + \cdots + y(a)$ pour tout $y \in F$.
Puisque $X_n = S(Y_n)$, la tribu $\mathcal{F}_n$ est incluse dans la tribu $\mathcal{G}_n$. On a donc,
pour toute fonction $f$ de $E$ dans $\mathbb{R}$,

$$
\mathrm{E}(f(X_{n+1}) | \mathcal{F}_n) = \mathrm{E}(\mathrm{E}(f(X_{n+1}) | \mathcal{G}_n) | \mathcal{F}_n).
$$

Par définition de X et puisque Y est une chaîne de Markov,

$$
\mathrm{E}(f(X_{n+1}) | \mathcal{G}_n) = \mathrm{E}(f(S(Y_{n+1})) | \mathcal{G}_n) = \mathrm{E}(f(S(Y_{n+1})) | Y_n).
$$

Enfin,

$$
\mathrm{E}(f(S(Y_{n+1})) | Y_n) = \left(1 - \frac{S(Y_n)}{a}\right) f(S(Y_n) + 1) + \left(\frac{S(Y_n)}{a}\right) f(S(Y_n) - 1).
$$

On a donc

$$
\begin{align*}
\mathrm{E}(f(X_{n+1}) | \mathcal{F}_n) &= \mathrm{E} \left( \left( 1 - \frac{X_n}{a} \right) f(X_n + 1) + \left( \frac{X_n}{a} \right) f(X_n - 1) \mid \mathcal{F}_n \right) \\
&= \left( 1 - \frac{X_n}{a} \right) f(X_n + 1) + \left( \frac{X_n}{a} \right) f(X_n - 1),
\end{align*}
$$

et donc $X$ est une chaîne de Markov de matrice de transition $\mathbf{P}$. Elle est irréductible puisque $i$ et $i+1$ communiquent pour tout $0 \le i \le a$ avec la convention $a+1=0$. Comme l'espace d'états est fini, la matrice est récurrente (positive). Pour tout $n$, $X_n$ et $X_{n+1}$ n'ont pas la même parité or $\mathbf{P}_{0,0}^2 > 0$ donc la chaîne est de période 2. Note : en général, l'image d'une chaîne de Markov par une fonction n'est pas une chaîne de Markov. $\square$

9.3 Pression au fil du temps

La pression au temps $n$ dans l'urne $A$ est de l'ordre de

$$
P_n := \frac{X_n}{a}.
$$

On s'intéresse ici aux deux premiers moments de cette variable aléatoire.
**Théorème 9.8** (Moyenne et variance de la pression). *Pour tout $n \in \mathbb{N}$, en notant $\alpha = 1 - 2/a$ et $\beta = 1 - 4/a$,*

$$
\begin{align*}
\mathrm{E}(P_n) &= \frac{1}{2} + \left( \mathrm{E}(P_0) - \frac{1}{2} \right) \alpha^n, \\
\mathrm{Var}(P_n) &= \frac{1 - \beta^n}{4a} + \mathrm{Var}(P_0)\beta^n - \left( \mathrm{E}(P_0) - \frac{1}{2} \right)^2 (\alpha^{2n} - \beta^n).
\end{align*}
$$

Démonstration. On raisonne par conditionnement :

$$
\begin{align*}
\mathrm{E}(X_{n+1} | X_n) &= (X_n + 1) \left(1 - \frac{X_n}{a}\right) + (X_n - 1) \frac{X_n}{a} \\
&= \alpha X_n + 1.
\end{align*}
$$

On en déduit, en prenant l'espérance, une relation de récurrence pour la suite
$(\mathbb{E}(P_n))_n$ qui donne immédiatement la formule pour la moyenne. De même,

$$
\begin{align*}
\mathrm{E}(X_{n+1}^2 | X_n) &= (X_n + 1)^2 \left(1 - \frac{X_n}{a}\right) + (X_n - 1)^2 \frac{X_n}{a} \\
&= \left(1 - \frac{4}{a}\right) X_n^2 + 2X_n + 1.
\end{align*}
$$

Ceci implique

$$
\begin{align*}
\mathrm{Var}(X_{n+1}) &= \beta\mathrm{Var}(X_n) + \frac{4}{a^2}\mathrm{E}(X_n)(a - \mathrm{E}(X_n)) \\
&= \beta\mathrm{Var}(X_n) + 1 - \frac{4}{a^2}\left(\mathrm{E}(X_0) - \frac{a}{2}\right)^2 \alpha^{2n}.
\end{align*}
$$

Un simple calcul fournit alors la formule pour la variance. $\square$

**Remarque 9.9** (Valeurs propres et vecteurs propres). *Soit $f$ la fonction sur $E$ définie par $f(x) = x - b$ où $b = a/2$, qui s'identifie avec le vecteur colonne $(-b, -b+1, \dots, b-1, b)^T$. Un calcul immédiat montre que $\mathbf{P}f = (1-2/a)f$. Ainsi, pour tout $n \ge 0$ et tout $x \in E$,*

$$
\mathrm{E}(X_n | X_0 = x) - b = \mathrm{E}(f(X_n) | X_0 = x) = \mathbf{P}^n f(x) = \alpha^n f(x) = \alpha^n (x - b).
$$

On retrouve la relation pour la moyenne du théorème 9.8 en prenant l'espérance par rapport à la loi de $X_0$. Plus généralement, on peut montrer que si $a = 2b$ est pair alors les valeurs propres de $\mathbf{P}$ sont les rationnels $(i/b)_{-b \le i \le b}$.

D'après le théorème 9.8, si le vide est fait dans l'urne A à l'instant initial,
alors $X_0 = 0$, $P_0 = 0$, et la pression moyenne $\mathrm{E}(P_n)$ au temps $n$ converge
géométriquement vers $1/2$ quand $n \to \infty$. D'autre part $\mathrm{Var}(P_n) \le 1/(4a)$
puisque $\mathrm{Var}(P_0) = 0$ et $\alpha^2 \ge \beta$, et l'inégalité de Tchebychev donne
$$P(|P_n - E(P_n)| \ge \varepsilon) \le \frac{1}{4a\varepsilon^2}.$$

Comme $a$ est très grand, la variable $P_n$ est quasiment déterministe, ce qui correspond à l'expérience. La figure 9.2 donne une trajectoire de $(P_n)_{n \ge 0}$ pour différentes conditions initiales.

<figure><img src="image_29.png" /><figcaption>Fig. 9.2. Évolution de la proportion de molécules dans l'urne A pour différentes valeurs initiales 0, 0.25, et 1, pour $a = 4000$ molécules.</figcaption></figure>

## 9.4 Fluctuations autour de la moyenne à l'équilibre

Supposons que le système soit à l'équilibre : $X$ suit la loi binomiale Bin($a, 1/2$). Alors, $P = X/a$ vérifie, en vertu du théorème limite central,

$$P\left(|P - 1/2| \ge \frac{x}{2\sqrt{a}}\right) \underset{a \text{ grand}}{\approx} 2\left(1 - \Phi(x)\right) \le \frac{2}{x\sqrt{2\pi}}e^{-x^2/2},$$

où $\Phi$ est la fonction de répartition de la loi gaussienne standard $N(0, 1)$.

On peut par ailleurs obtenir une estimation non asymptotique de cette probabilité en utilisant l'inégalité de Chernoff. Plus précisément, si $(Z_i)_{1 \le i \le n}$ sont des variables aléatoires indépendantes et identiquement distribuées de loi de Rademacher $(1/2)(\delta_{-1} + \delta_1)$ alors, pour tout $r \ge 0$ et tout $\lambda > 0$,
$$ \mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n} Z_i \ge r\right) = \mathbb{P}\left(\exp\left(\lambda \sum_{i=1}^{n} Z_i\right) \ge e^{\lambda nr}\right) \le e^{-\lambda nr} \mathrm{ch}(\lambda)^n. $$

Comme $\mathrm{ch}(\lambda) \le \exp(\lambda^2/2)$ (découle par exemple d'un développement en série), on obtient, après une optimisation en $\lambda$ qui conduit à choisir $\lambda = r$,

$$ \mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n} Z_i \ge r\right) \le e^{-nr^2/2}. $$

Ceci assure que

$$ \mathbb{P}(|P - 1/2| \ge r) \le 2e^{-ar^2/2}. $$

Si $a = 6 \times 10^{23}$ alors, en prenant $r = 10^{-9}$, la probabilité de trouver une fluc- tuation supérieure au milliardième est inférieure à $2e^{-3 \times 10^5}$. Cette probabilité extrêmement petite rend l'événement inaccessible.

## 9.5 Quand reverra-t-on le vide ?

La chaîne d'Ehrenfest $(X_n)_{n \ge 0}$ est irréductible à espace d'états fini, et elle visite donc tous les états, et le temps d'atteinte de chaque état est inté- grable (en particulier fini p.s.). Ceci peut sembler paradoxal si l'on pense à l'irréversibilité macroscopique. Essayons donc d'estimer ces temps d'atteinte moyens. Le plus simple est de s'intéresser tout d'abord aux temps de retour. Sur l'événement $\{X_0 = i\}$, le temps de retour $T_{ii}$ en $i$ est la variable aléatoire

$$ T_{ii} := \inf\{n \ge 1 : X_n = i\}. $$

**Théorème 9.10 (Temps de retours moyens).** Si $a$ est pair et $b := a/2$ alors

$$ E(T_{00} | X_0 = 0) = 2^a \quad \text{et} \quad E(T_{bb} | X_0 = b) \sim \sqrt{\pi b}. $$

*Démonstration.* Notons que le premier résultat découle du lemme 9.1. On sait que $E(T_{ii} | X_0 = i)$ est égale à $1/\pi(i)$ où $\pi$ est la loi invariante, qui n'est rien d'autre que $\mathrm{Bin}(a, 1/2)$. On a $\pi(i) = 2^{-a} \binom{a}{i}$, et le résultat est alors immédiat pour $i = 0$ et découle de la formule de Stirling pour $i = b$. $\square$

Il faut donc attendre un temps immense avant que l'urne A initialement vide, ne le redevienne ! On souhaite à présent plus de précision et obtenir des estimations des temps d'atteinte de 0 et b partant de 0 ou b. Soit

$$ m_{i,j} := \mathbb{E}(\inf\{n \ge 0 : X_n = j\} | X_0 = i) $$

le temps moyen d'atteinte de j par la chaîne partant de i. Partant de $i+1$, la chaîne doit passer par i pour aller en 0, donc pour $0 \le i < a$,

$$ m_{i+1,0} = m_{i+1,i} + m_{i,0}. $$
Par symétrie, pour $0 \le i < a$,

$$m_{i,i+1} = m_{a-i,a-i-1} = m_{a-i,0} - m_{a-i-1,0}.$$

Ainsi les réels $(m_{i,j})_{ij}$ sont déterminés par les réels $(m_{i,0})_i$. Remarquons à présent que $X_n = 0$ si et seulement si $Y_n = (0, \dots, 0)$. Donc le temps moyen que met $X$ à se rendre de $i$ à 0 est le temps moyen que met $Y$ pour aller de n'importe quel élément de $F$ situé à une distance $i$ de $(0, \dots, 0)$ à $(0, \dots, 0)$, autrement dit $m_i$. Le théorème 9.2 assure alors que pour $0 \le i \le a - 1$,

$$m_{i,i+1} = m_{a-i,0} - m_{a-i-1,0} = \sum_{k=1}^{a-i} \mathbf{Q}_{a-k}^{a} - \sum_{k=1}^{a-i-1} \mathbf{Q}_{a-k}^{a} = \mathbf{Q}_{i}^{a}.$$

**Corollaire 9.11 (Temps d'atteinte moyens).** On a

$$m_{i,j} = \begin{cases} 2^a / \binom{a}{i} & \text{si } i = j, \\ \displaystyle\sum_{k=i}^{j-1} \mathbf{Q}_k^a & \text{si } i < j, \\ \displaystyle\sum_{k=a-i}^{a-j-1} \mathbf{Q}_k^a & \text{si } i > j. \end{cases}$$

**Théorème 9.12 (Temps d'atteinte moyens).** Si $a$ est pair et $b := a/2$ alors

$$a\left(2\log(2) - \frac{1}{2}\right) - 2 \le m_{0,b} \le a\left(\frac{1}{4} + \frac{1}{2}\log(a)\right),$$

et

$$2^a \left(1 + \frac{1}{a-1}\right) \le m_{0,a} \le 2^a \left(1 + \frac{2}{a-1}\right).$$

*Preuve du second encadrement du théorème 9.12.* Pour tout $0 \le i \le a - 1$,

$$m_{i,i+1} + m_{i+1,i} = \frac{1}{\binom{a-1}{i}} \left( \sum_{k=0}^{i} \binom{a}{k} + \sum_{k=i+1}^{a} \binom{a}{k} \right) = \frac{2^a}{\binom{a-1}{i}}.$$

Donc

$$m_{0,b} + m_{b,0} = 2^a \sum_{i=0}^{b-1} \frac{1}{\binom{a-1}{i}}.$$

De plus, puisque $m_{b,0} = m_{b,a}$, on a

$$m_{0,a} = m_{0,b} + m_{b,a} = 2^a \sum_{i=0}^{b-1} \frac{1}{\binom{a-1}{i}}.$$

En séparant les deux premiers termes de la somme des autres, on a
$$1 + \frac{1}{a-1} \le \sum_{i=0}^{b-1} \frac{1}{\binom{a-1}{i}} \le 1 + \frac{1}{a-1} + \frac{b-1}{\binom{a-1}{2}}$$

ce qui fournit l'encadrement souhaité. □

Le théorème 9.12 assure que le temps mis par le système pour passer d'un état de total déséquilibre à celui d'équilibre parfait est extrêmement négligeable devant le temps mis pour l'évolution inverse. Pour un nombre de molécules égal à 100, $m_{0,50}$ est majoré par 256 tandis que $m_{50,0}$ est de l'ordre de $10^{30}$ soit quelques mille milliards de milliards de milliards.

## 9.6 Pour aller plus loin

Le temps de mélange pour la marche aléatoire sur le cube discret vérifie une propriété de convergence abrupte autour de $(a/2) \log(a)$ de fenêtre $a$: dans un intervalle de longueur de l'ordre de $a$ autour de $(a/2) \log a$ la distance en variation totale entre $\mathbf{R}^n(x, \cdot)$ et $\pi$ passe d'un voisinage de 1 à un voisinage de 0. La borne supérieure est fournie par le théorème 9.4, qu'on peut compléter par l'estimation suivante qui se trouve dans le livre David Levin, Yuval Peres, et Elizabeth Wilmer [LPW09] : si $n = (1/2)a \log(a) - ca$ avec $c > 0$, alors

$$\max_{x \in F} d_{\text{VT}}(\mathbf{R}^n(x, \cdot), \pi) \ge 1 - 8e^{-(2c+1)}.$$

Ce résultat semble provenir d'un travail de Persi Diaconis et Mehrdad Shahshahi [DS87]. Les valeurs propres et vecteurs propres de la matrice de transition $\mathbf{P}$ de la chaîne d'Ehrenfest sont obtenus par Mark Kac dans [Kac47]. On peut également les obtenir à partir de ceux de la matrice de transition $\mathbf{R}$ du modèle microscopique (marche sur le cube discret) qui sont décrits dans le livre de Levin, Peres, et Wilmer [LPW09, exemple 12.15].

La relation classique entre espérance du temps de retour et loi invariante utilisée dans la preuve du théorème 9.10 est établie par exemple dans le livre de James Norris [Nor98a]. Les théorèmes 9.2 et 9.12 sont tirés du livre de John Kemeny et Laurie Snell [KS76].

Les molécules passent très vite d'une urne à l'autre. Ainsi, le temps microscopique est différent du temps macroscopique. Considérons qu'une unité de temps macroscopique corresponde à $a$ changements d'urnes. Quelle est alors l'évolution de la pression dans cette nouvelle échelle de temps ? Le théorème 9.8 entraîne immédiatement que, pour tout $t > 0$,

$$E(P_{[at]}) \xrightarrow{a \to \infty} \frac{1}{2} + \left(E(P_0) - \frac{1}{2}\right) e^{-2t} \quad \text{et} \quad E(V_{[at]}) \xrightarrow{a \to \infty} \text{Var}(X_0)e^{-4t},$$

puisque $[at]/a$ tend vers $t$. Cette convergence est la partie émergée d'un iceberg : convenablement renormalisée en temps et en espace, la chaîne d'Ehrenfest converge en loi vers le processus de diffusion d'Ornstein-Uhlenbeck, comme présenté dans le chapitre 27.
<figure><img src="image_30.png" /></figure>
# Records, extrêmes, et recrutements

**Mots-clés.** Extrême; statistique d'ordre.

**Outils.** Processus de Markov déterministe par morceaux; martingale; permutation aléatoire; théorème des séries centrées; théorème de Lindeberg-Lévy; loi exponentielle; loi de Pareto; transformée de Laplace.

**Difficulté.** **

Un cabinet de chasseurs de têtes souhaite constituer une équipe de choc grâce à un recrutement au fil de l'eau : chaque candidat obtient une note après son entretien et les recruteurs décident dans l'instant de l'engager ou pas. Le cahier des charges est de recruter une équipe à la fois excellente et fournie. Dans ce chapitre, on compare deux politiques de recrutement différentes : une personne est recrutée si elle obtient une note supérieure à toutes celles des candidats précédents pour la première et à la moyenne de ces notes pour la seconde. La deuxième procédure, plus souple, devrait conduire à un recrutement plus important mais peut-être de qualité moindre. On modélise les notes des candidats successifs par une suite $(X_n)_{n \ge 1}$ de variables aléatoires réelles indépendantes et de même loi $\mu$ de fonction de répartition $F$. Afin d'éviter les cas particuliers, nous supposerons que $F$ est continue sur $\mathbb{R}$. L'extrémité droite du support de $\mu$ est notée $x_F$ :

$$x_F = \sup \{x \in \mathbb{R} : F(x) < 1\} \in \mathbb{R} \cup \{+\infty\}.$$

## 10.1 Élitisme

Pour tout $n \ge 1$, le rang relatif de $X_n$ dans $X_1, \dots, X_n$ est donné par

$$R_n = 1 + \sum_{k=1}^{n-1} \mathbb{1}_{\{X_n > X_k\}}.$$
Le candidat $n$ est recruté si son rang vaut $n$, c'est-à-dire si sa note bat le précédent record. Le nombre de records, ou de personnes recrutées, jusqu'au temps $n$ est donné par

$$Z_n = \sum_{i=1}^{n} 1_{\{R_i=i\}}.$$

On définit également les instants de record par récurrence :

$$T_1 = 1 \text{ et, pour } j \ge 1, \quad T_{j+1} = \inf \{n > T_j : X_n > X_i \text{ pour } i < n\}.$$

**Lemme 10.1 (Lois des rangs).** *Les variables aléatoires $(R_n)_{n \ge 1}$ sont indépendantes et, pour tout $n \in \mathbb{N}^*$, $R_n$ suit la loi uniforme sur $\{1, \dots, n\}$,*

$$E(Z_n) = \sum_{k=1}^{n} \frac{1}{k} = \log(n) + \gamma + \mathcal{O}(1/n),$$

et

$$\mathrm{Var}(Z_n) = \sum_{k=1}^{n} \frac{k-1}{k^2} = \log(n) + \gamma - \frac{\pi^2}{6} + \mathcal{O}(1/n),$$

où $\gamma$ désigne la constante d'Euler.

*Démonstration.* La fonction de répartition $F$ de $\mu$ est continue donc pour tous entiers $i$ et $j$ distincts, $P(X_i \neq X_j) = 1$. En particulier, pour $n \ge 1$, presque sûrement, il existe une unique permutation (aléatoire) $\sigma$ de $\{1, \dots, n\}$ telle que

$$X_{\sigma(1)} < \dots < X_{\sigma(n)}.$$

De plus la loi de cette permutation est la loi uniforme sur $S_n$. Enfin, à un vecteur $(R_1, \dots, R_n)$ on associe la permutation $\sigma$ de manière bijective, d'où

$$P(R_1 = r_1, \dots, R_n = r_n) = \frac{1}{n!}$$

pour tout $(r_1, \dots, r_n)$ vérifiant $1 \le r_i \le i$ pour $1 \le i \le n$. La v.a. $Z_n$ est donc la somme de $n$ v.a. de loi de Bernoulli indépendantes de paramètres 1, 1/2, ..., 1/n. Les expressions de son espérance et de sa variance s'en déduisent. $\square$

Notons que la suite $(Z_n)_{n \ge 1}$ a la même loi que la suite $(|B_n|)_{n \ge 1}$ du processus des restaurants chinois du chapitre 14, voir également le théorème 13.5 sur la suite $(K_n)_{n \ge 1}$ du nombre d'allèles dans un échantillon obtenue lors de l'étude de la généalogie du modèle de Wright-Fisher.

**Théorème 10.2 (Comportement asymptotique du nombre de records).**

$$\frac{Z_n}{\log(n)} \xrightarrow{n \to \infty} 1 \quad \text{et} \quad \frac{Z_n - \log(n)}{\sqrt{\log(n)}} \xrightarrow{n \to \infty} \mathcal{N}(0, 1).$$
**Démonstration.** La convergence presque sûre repose sur le théorème des séries centrées que nous rappelons ci-dessous.

**Théorème 10.3 (Séries centrées).** Soit $(Q_n)_{n \ge 1}$ une suite de variables aléatoires indépendantes et centrées et une suite (déterministe) $(b_n)_{n \ge 1}$ strictement croissante non bornée de réels strictement positifs.

$$ Si \sum_{n=1}^{\infty} \frac{\mathbb{E}(Q_n^2)}{b_n^2} < \infty, \quad alors \quad \frac{1}{b_n} \sum_{k=1}^{n} Q_k \xrightarrow{n \to \infty} 0. $$

Reste à présent à appliquer ce théorème à la suite $(Q_n)_{n \ge 1}$ définie par $Q_n = \mathbb{1}_{\{R_n=n\}} - 1/n$ pour $n \ge 1$ avec $b_n = \log(n)$ en remarquant

$$ \frac{\mathrm{Var}(Q_n)}{b_n^2} = \frac{n-1}{n^2 \log(n)^2} $$

est le terme général d'une série convergente.

La deuxième partie du théorème découle du théorème limite central de Lindeberg-Lévy rappelé ci-dessous. Il s'agit d'une extension du théorème limite central aux suites de variables aléatoires indépendantes de carré intégrable, pas forcément de même loi. La condition de Lindeberg signifie que la variance totale est asymptotiquement bien répartie.

**Théorème 10.4 (Lindeberg-Lévy).** Soit $(Y_i)_{i \ge 1}$ une suite de v.a.r. indépendantes de carré intégrable et centrées. Notons $S_n = Y_1 + \dots + Y_n$ et $s_n$ son écart-type. Si la condition de Lindeberg est vérifiée : pour tout $\varepsilon > 0$,

$$ \lim_{n \to \infty} \frac{1}{s_n^2} \sum_{i=1}^{n} \mathbb{E}[Y_i^2 \mathbb{1}_{\{|Y_i| > \varepsilon s_n\}}] = 0, $$

alors

$$ \frac{S_n}{s_n} \xrightarrow{n \to \infty} \mathcal{N}(0, 1). $$

On applique ce théorème aux variables $Y_n = \mathbb{1}_{\{R_n=n\}} - 1/n$ pour $n \ge 1$ en remarquant que $\mathbb{1}_{\{|Y_i| > \varepsilon s_n\}} = 0$ pour $n$ assez grand. $\square$

**Remarque 10.5 (Avec les martingales).** Le théorème 10.2 peut également être démontré grâce à la loi des grands nombres pour les martingales et au théorème limite central pour les martingales. Posons, pour $n \ge 1$,

$$ M_n = Z_n - \sum_{k=1}^{n} \frac{1}{k} = \sum_{k=1}^{n} \left( \mathbb{1}_{\{R_k=k\}} - \frac{1}{k} \right). $$

Alors $(M_n)_{n \ge 1}$ est une martingale de carré intégrable de processus croissant

$$ \langle M \rangle_n = \sum_{k=1}^{n} \frac{k-1}{k^2} = \mathcal{O}(\log(n)). $$

Ainsi, $M_n = o(\log(n))$ presque sûrement. La convergence en loi se déduit également du théorème limite central pour $(M_n)_{n \ge 1}$.
## 10.2 Au-dessus de la moyenne

On construit alors par récurrence les suites $(Y_n)_{n \ge 1}$, $(\bar{Y}_n)_{n \ge 1}$ et $(T_n)_{n \ge 1}$ :

$$Y_1 = X_1, \quad \bar{Y}_1 = Y_1 \quad \text{et} \quad T_1 = 1,$$

puis, pour tout $n \ge 1$, par

$$T_{n+1} = \inf \{i > T_n : X_i > \bar{Y}_n\}, \quad Y_{n+1} = X_{T_{n+1}} \quad \text{et} \quad \bar{Y}_{n+1} = \frac{1}{n+1} \sum_{i=1}^{n+1} Y_i.$$

Pour tout $n \ge 1$, $Y_n$ est la performance du $n$-ième candidat retenu, $T_n$ est le nombre de candidats auditionnés pour en retenir $n$ et $\bar{Y}_n$ est la moyenne des performances des $n$ premiers candidats retenus.

**Théorème 10.6** (Vers la qualité maximale).

$$\bar{Y}_n \xrightarrow{n \to \infty} x_F.$$

*Démonstration.* La suite $(\bar{Y}_n)_{n \ge 1}$ est croissante par construction. Elle converge donc vers une variable aléatoire $Y$ à valeurs dans $\bar{\mathbb{R}}$. Nous allons utiliser la formule suivante : $\bar{Y}_k = \sum_{j=1}^k (Y_j - \bar{Y}_{j-1})/j$. Soit à présent $x < x_F$ et $\varepsilon > 0$ tel que $x + \varepsilon < x_F$. Si $x_F = +\infty$, on peut prendre $\varepsilon = 1$. Sur l'évènement $\{Y \le x\}$ on a, puisque $\bar{Y}_j \le Y \le x$,

$$Y = \sum_{j=1}^{\infty} \frac{Y_j - \bar{Y}_{j-1}}{j} \ge \sum_{j=1}^{\infty} \frac{(Y_j - x)^+}{j} \ge \sum_{j=1}^{\infty} \frac{\varepsilon \mathbf{1}_{\{Y_j > x+\varepsilon\}}}{j}.$$

Posons, pour $k \ge 1$, $M_k = \#\{j \le k : Y_j > x + \varepsilon\}$ et $N_k = \#\{j \le k : X_j > x + \varepsilon\}$. Alors $M_k \ge N_k$ pour tout $k \ge 1$ et $\mathbf{1}_{\{Y_j > x+\varepsilon\}} = M_j - M_{j-1}$. On a donc

$$\sum_{j=1}^{\infty} \frac{\mathbf{1}_{\{Y_j > x+\varepsilon\}}}{j} = \sum_{j=1}^{\infty} \frac{M_j - M_{j-1}}{j} = \sum_{j=1}^{\infty} \frac{M_j}{j(j+1)} \\ \ge \sum_{j=1}^{\infty} \frac{N_j}{j(j+1)} \ge \sum_{j=1}^{\infty} \frac{N_j - N_{j-1}}{j} = \sum_{j=1}^{\infty} \frac{\mathbf{1}_{\{X_j > x+\varepsilon\}}}{j}.$$

Il reste à remarquer que presque sûrement la série de terme général $\mathbf{1}_{\{X_j > x+\varepsilon\}}/j$ diverge puisque la série centrée associée converge mais pas la série des espérances. Elle diverge notamment sur l'ensemble $\{Y \le x\}$ ce qui implique que cet ensemble est de probabilité nulle. $\square$

On s'intéresse alors au nombre de candidats qu'il faudra auditionner pour former une équipe de taille donnée. Le résultat suivant permet de conclure sous une hypothèse assez souvent vérifiée en pratique sur la convergence de la suite $(\bar{Y}_n)_{n \ge 1}$.
**Théorème 10.7** (Comportement asymptotique des auditionnés). Soit la suite de variables aléatoires $(P_n)_{n \ge 1}$ définie par $P_1 = 1$ et, pour $n \ge 2$, $P_n = 1 - F(\bar{Y}_{n-1})$. S'il existe un réel $\alpha > 0$ et une variable aléatoire $W$ tels que $P(0 < W < \infty) = 1$ et $n^\alpha P_n \xrightarrow{n \to \infty} W$ alors

$$ \frac{T_n}{n^{\alpha+1}} \xrightarrow{n \to \infty} \frac{1}{(\alpha+1)W}. $$

*Démonstration.* L'idée est d'appliquer le théorème 10.3 conditionnellement à la tribu $\bar{\mathcal{F}}$ engendrée par la suite $(\bar{Y}_n)_{n \ge 1}$ en posant, pour $n \ge 1$,

$$ b_n = \sum_{i=1}^{n} \frac{1}{P_i} \quad \text{et} \quad Q_n = T_n - T_{n-1} - \frac{1}{P_n} \quad \text{avec} \quad T_0 = 0. $$

La suite $(b_n)_{n \ge 1}$ est strictement croissante et tend vers $+\infty$ puisque $1/P_i$ est de l'ordre de $i^\alpha/W$. Plus précisément, une comparaison série-intégrale donne

$$ b_n \sim \frac{n^{\alpha+1}}{(\alpha+1)W}. $$

D'autre part, $T_n - T_{n-1}$ suit la loi géométrique de paramètre $P_n$ donc

$$ \mathbb{E}\left(\frac{Q_n^2}{b_n^2} \mid \bar{\mathcal{F}}\right) = \frac{1-P_n}{b_n^2 P_n^2} \underset{n \to \infty}{\sim} \frac{(\alpha+1)^2}{n^2} $$

qui est le terme général d'une série convergente. Le théorème des séries centrées (théorème 10.3) assure donc que

$$ \frac{1}{b_n} \sum_{k=1}^{n} Q_k = \frac{T_n}{b_n} - 1 \xrightarrow{n \to \infty} 0. $$

L'équivalent de $b_n$ donne le résultat annoncé. $\square$

## 10.3 Cas de la loi exponentielle

**Théorème 10.8** (Évolution de la moyenne et de l'effectif). Supposons que les variables aléatoires $(X_n)_{n \ge 1}$ soient distribuées selon la loi exponentielle de paramètre $\lambda$. Alors il existe une variable aléatoire $G$ de loi de Gumbel de fonction de répartition $x \mapsto \exp(-e^{-x})$ telle que

$$ \lambda \bar{Y}_n - \log(n) \xrightarrow{n \to \infty} G, \quad \lambda \bar{Y}_n - \log(n) \xrightarrow{n \to \infty} G, \quad \text{et} \quad \frac{T_n}{n^2} \xrightarrow{n \to \infty} \frac{e^G}{2}. $$
<figure><img src="image_31.png" /><figcaption>Fig. 10.1. Convergence p.s. de la suite $(\overline{Y}_n - \log(n))_{n \ge 2}$.</figcaption></figure>

*Démonstration.* Quitte à considérer la suite $(\lambda Y_n)_{n \ge 1}$, on peut supposer que $\lambda = 1$. Avec la convention $\overline{Y}_0 = 0$, posons, pour tout $n \ge 0$,

$$\eta_n = Y_{n+1} - \overline{Y}_n = (n+1)(\overline{Y}_{n+1} - \overline{Y}_n).$$

On a donc $\overline{Y}_n = \sum_{j=1}^n \eta_{j-1}/j$. Notons $\mathcal{F}_n$ la tribu engendrée par $(Y_1, \dots, Y_n)$. Pour tout $n \ge 1$, la loi conditionnelle de $\eta_n$ sachant $\mathcal{F}_n$ est égale à la loi de $Y_{n+1} - \overline{Y}_n$ sachant $\overline{Y}_n$ et $\{Y_{n+1} > \overline{Y}_n\}$, d'où, en utilisant la propriété d'absence de mémoire des lois exponentielles,

$$P(\eta_n > x | \mathcal{F}_n) = P(Y_{n+1} > \overline{Y}_n + x | \overline{Y}_n, Y_{n+1} > \overline{Y}_n) = e^{-x}.$$

Les variables aléatoires $(\eta_n)_{n \ge 0}$ sont donc indépendantes et de même loi exponentielle de paramètre 1. Avec $h_n = \sum_{k=1}^n \frac{1}{k}$, on a

$$\overline{Y}_n - h_n = \sum_{k=1}^{n} \frac{\eta_{k-1} - 1}{k}.$$

La suite $(\overline{Y}_n - h_n)_{n \ge 1}$ est une martingale (sommes partielles de variables aléatoires indépendantes intégrables et centrées). Elle est bornée dans $L^2$ car
<figure><img src="image_32.png" /><figcaption>Fig. 10.2. Histogramme de la limite presque sûre de la suite aléatoire $(\bar{Y}_n - \log n)_n$ et densité de la loi de Gumbel.</figcaption></figure>

$$ \mathrm{E}((\bar{Y}_n - h_n)^2) = \sum_{k=1}^{n} \frac{\mathrm{E}((\eta_{k-1} - 1)^2)}{k^2} \le \mathrm{E}((\eta_0 - 1)^2) \frac{\pi^2}{6} < \infty. $$

Il existe donc une variable aléatoire $\tilde{G}$ telle que $\bar{Y}_n - h_n$ converge vers $\tilde{G}$ dans $L^2$ et presque sûrement. Nous allons à présent établir que $(\bar{Y}_n - \log(n))_{n \ge 1}$ converge en loi vers la loi de Gumbel. Il en découlera que $\tilde{G} = G - \gamma$ où G suit la loi de Gumbel et $\gamma = \lim_{n \to \infty} (h_n - \log(n))$ est la constante d'Euler.

La transformée de Laplace de $\bar{Y}_n$ notée $H_n$ vaut, pour tout $\theta \ge 0$,

$$ H_n(\theta) = \mathrm{E}(e^{-\theta \bar{Y}_n}) = \prod_{k=1}^{n} \frac{k}{k+\theta} = n! \left( \prod_{k=1}^{n} (k+\theta) \right)^{-1}. $$

Soit à présent $(Z_n)_{n \ge 1}$ une suite de v.a.r. i.i.d. de loi exponentielle de paramètre 1. Pour tout $n \ge 1$, soit $Z_{(n)} = \max(Z_1, \dots, Z_n)$. Calculons la transformée de Laplace $L_n$ de $Z_n$. Or $Z_{(n)}$ a pour densité $ne^{-z}(1-e^{-z})^{n-1}1_{[0, +\infty)}(z)$ puisque $\mathbb{P}(Z_{(n)} \le z) = \mathbb{P}(Z_1 \le z) \cdots \mathbb{P}(Z_n \le z) = (1-e^{-z})^n 1_{[0, +\infty)}(z)$. Grâce à une intégration par parties, on a

$$ L_n(\theta) := \mathrm{E}(e^{-\theta Z_{(n)}}) $$
$$
\begin{align*}
&= \int_{0}^{\infty} e^{-\theta z} n e^{-z} (1 - e^{-z})^{n-1} dz \\
&= \frac{n}{1+\theta} \int_{0}^{\infty} e^{-\theta z} (n-1) e^{-2z} (1 - e^{-z})^{n-2} dz \\
&= \frac{n}{1+\theta} L_{n-1}(\theta+1),
\end{align*}
$$

et une récurrence immédiate donne

$$
L_n(\theta) = \prod_{k=1}^{n} \frac{n+1-k}{k+\theta} = n! \left( \prod_{k=1}^{n} (k+\theta) \right)^{-1} .
$$

On en déduit que $\bar{Y}_n$ et $Z_{(n)}$ ont même loi. En fait cette propriété ne constitue qu'une partie du lemme de Rényi 11.3. Ainsi, pour tout réel $x$,

$$
\mathbb{P}(\bar{Y}_n - \log(n) \le x) = \mathbb{P}(Z_{(n)} - \log(n) \le x) = \left(1 - \frac{e^{-x}}{n}\right)^n \underset{n \to \infty}{\longrightarrow} e^{-e^{-x}},
$$

ce qui conclut la preuve des deux premiers points. Pour le dernier point, on a

$$
n P_n = n(1 - F(\bar{Y}_{n-1})) = e^{-(\bar{Y}_n - \log(n))} \xrightarrow{n \to \infty} e^{-G},
$$

ce qui donne $T_n/n^2 \xrightarrow{n \to \infty} e^G/2$ par le théorème 10.7. $\square$

10.4 Cas de la loi de Pareto

**Théorème 10.9** (Évolution de la moyenne). *Supposons que les variables aléatoires* $(X_n)_{n \ge 1}$ *soient distribuées selon la loi de Pareto de paramètre* $(\theta, \beta)$ *avec* $\theta > 0$ *et* $\beta > 1$, *c'est-à-dire qu'elles possèdent la fonction de répartition suivante :*

$$
F(x) = \left[ 1 - \left( \frac{\theta}{x} \right)^{\beta} \right] 1_{[\theta, +\infty[}(x).
$$

Alors, il existe une variable aléatoire Y telle que $\mathbb{P}(0 < Y < \infty) = 1$ et

$$
\frac{\bar{Y}_n}{n^{1/(\beta-1)}} \xrightarrow{n \to \infty} Y.
$$

Démonstration. Si $X$ suit la loi de Pareto de paramètre $(\theta, \beta)$ alors $X/\theta$ suit la loi de Pareto de paramètre $(1, \beta)$. On supposera donc dans la suite de la preuve que $\theta = 1$. De plus, si $X$ suit la loi de Pareto de paramètre $(1, \beta)$, alors la loi de $X$ sachant que $X > c$ est la loi de $cX$. En notant $\mathcal{F}_n$ la tribu engendrée par $(Y_1, \dots, Y_n)$, on a, pour tout $x > 1$,

$$
\mathbb{P}(Y_n > x \bar{Y}_{n-1} | \mathcal{F}_{n-1}) = \mathbb{P}(X > x \bar{Y}_{n-1} | X > \bar{Y}_{n-1}, \bar{Y}_{n-1}) = x^{-\beta}.
$$
où X est une variable aléatoire de loi de Pareto de paramètre $\beta$ indépendante de $\mathcal{F}_{n-1}$. En d'autres termes, la loi de $Y_n/\bar{Y}_{n-1}$ sachant $\mathcal{F}_{n-1}$ est la loi de Pareto de paramètre $\beta$. En conséquence,

$$ \mathbb{E}(\bar{Y}_n | \mathcal{F}_{n-1}) = \frac{(n-1)\bar{Y}_{n-1} + \mathbb{E}(Y_n | \bar{Y}_{n-1})}{n} = \frac{n-1 + \mathbb{E}(X)}{n} \bar{Y}_{n-1}. $$

Posons, pour $n \ge 1$,

$$ a_n = \frac{n-1 + \mathbb{E}(X)}{n} = 1 + \frac{1}{n(\beta-1)} \quad \text{et} \quad b_n = \frac{1}{\prod_{j=1}^{n} a_j}. $$

La suite $(M_n)_{n \ge 1}$ définie par $M_n = b_n \bar{Y}_n$ est une martingale positive d'espérance 1. Elle converge donc presque sûrement vers une variable aléatoire $M$ positive de $L^1$. Nous omettons ici la preuve assez technique du fait que $\mathbb{P}(0 < M < \infty) = 1$. Enfin, on remarque que

$$ b_n = \exp\left(-\sum_{j=1}^{n} \log\left(1 + \frac{1}{(\beta-1)j}\right)\right) = \exp\left(-\frac{1}{\beta-1} \sum_{j=1}^{n} \frac{1}{j} + o(1)\right). $$

En d'autres termes, $n^{1/(\beta-1)}b_n$ converge vers un réel strictement positif. $\square$

**Corollaire 10.10 (Temps de recrutement).** *Il existe une variable aléatoire T à valeurs dans $(0, \infty)$ telle que*

$$ \frac{T_n}{n^{(2\beta-1)/(\beta-1)}} \xrightarrow{n \to \infty} T. $$

*Démonstration.* Pour tout $n \ge 1$, on a $P_n = (\bar{Y}_{n-1})^{-\beta}$. Le résultat précédent assure que $n^{\beta/(\beta-1)}P_n$ converge presque sûrement vers une variable aléatoire $W$ vérifiant $\mathbb{P}(0 < W < \infty) = 1$. On conclut alors avec le théorème 10.7. $\square$

## 10.5 Pour aller plus loin

On peut aussi étudier le nombre $M_n$ de candidats retenus après $n$ auditions et la moyenne $A_n$ des performances des candidats retenus après $n$ auditions.

Ce chapitre est essentiellement basé sur un article de John Preater [Pre00] pour le cas de la loi exponentielle, et sur un article de Abba Krieger, Moshe Pollak, et Ester Samuel-Cahn [KPSC08] pour le cas général. Les preuves du théorème des séries centrées (théorème 10.3) et du théorème de Lindeberg-Lévy (théorème 10.4) se trouvent par exemple dans l'incontournable livre de William Feller [Fel68].

On peut compléter le théorème 10.2 avec des inégalités de déviation non asymptotiques, à la différence du théorème de convergence en loi. Comme les variables aléatoires $(R_n)_{n \ge 1}$ sont à valeurs dans $[0, 1]$, l'inégalité de (Wassily) Hoeffding assure que, pour tout $r > 0$,
$$P\left(Z_n - \sum_{k=1}^{n} \frac{1}{k} \geq r\right) \leq \exp\left(-\frac{2r^2}{n}\right).$$

Cette inégalité n'est en l'occurrence pas très précise car elle ne prend pas bien en compte le fait que les variances des incrémentes de $(Z_n)_{n \ge 1}$ tendent vers 0. On peut obtenir un résultat plus fin pour la déviation inférieure grâce à l'inégalité de (Andreas) Maurer pour les sommes de variables aléatoires positives [Mau03], qui donne, pour tout $r > 0$,

$$P\left(Z_n - \sum_{k=1}^{n} \frac{1}{k} \leq -r\right) \leq \exp\left(-\frac{r^2}{2\log(n)}\right).$$

Citons enfin le célèbre problème des secrétaires<fnref n="1" /> dont on trouve une description savoureuse dans un article de Thomas Ferguson [Fer89]. Le contexte est le même que dans notre chapitre mais on ne recrute qu'une personne parmi $n$ candidats. La question est de déterminer la stratégie fournissant avec la probabilité maximale le meilleur candidat. La première étape est de montrer que la stratégie optimale est nécessairement de la forme suivante : on fixe $r$ entre 1 et $n$, on recale les $r-1$ premiers candidats puis on sélectionne le premier candidat qui obtient une meilleure note que ses prédécesseurs. Avec une telle stratégie, on recrute le meilleur candidat au rang $j$ (entre $r$ et $n$) si et seulement si le meilleur candidat arrive en position $j$ et si les candidats numérotés de $r$ à $j-1$ sont moins bons que l'un des $r-1$ premiers. Ainsi,

$$P(j \text{ est le meilleur et il est sélectionné}) = \frac{1}{n} \frac{r-1}{j-1}$$

et la probabilité de recruter le meilleur candidat avec cette stratégie est

$$p_n(r) = \sum_{j=r}^{n} \frac{1}{n} \frac{r-1}{j-1} = \frac{r-1}{n} \sum_{j=r}^{n} \frac{1}{j-1}.$$

Pour $x \in ]0, 1[$,

$$p_n(\lfloor xn \rfloor) = \frac{\lfloor xn \rfloor}{n} \sum_{j=\lfloor xn \rfloor}^{n} \frac{1}{j-1} \underset{n \to \infty}{\longrightarrow} -x \log(x).$$

La probabilité maximale vaut donc $1/e$ pour un choix optimal de $r \approx n/e$.

<footnote n="1">1. Ce problème d'arrêt optimal est un grand classique de la théorie du contrôle stochastique et de la théorie de la décision. Il est également connu sous le nom du problème du mariage (entre autres).</footnote>
<figure><img src="image_33.png" /><figcaption>**Fig. 10.3.** Convergence de la probabilité optimale du meilleur choix (croix) et de la proportion idéale (trait continu) vers 1/e.</figcaption></figure>
<figure><img src="image_34.png" /></figure>
File d'attente M/M/Infini

**Mots-clés.** Loi exponentielle; loi binomiale; loi de Poisson; file d'attente; loi invariante.

**Outils.** Chaîne de Markov à temps continu; processus de Poisson; couplage; distance de Wasserstein; semi-groupe; générateur infinitésimal; fonction génératrice; inégalité de Markov.

**Difficulté.** **

Les files d'attente<fnref n="1" /> font partie des modèles aléatoires les plus répandus et les plus utiles. Le cas le plus simple à décrire est sans doute le suivant : des clients font la queue devant un guichet appelé serveur. Les durées qui séparent les arrivées des clients successifs sont modélisées par des v.a.r. i.i.d. de loi exponentielle de paramètre $\lambda$, tandis que les durées de traitement des clients successifs par le serveur sont modélisées par des v.a.r. i.i.d. de loi exponentielle de paramètre $\mu$. Le choix de la loi exponentielle est justifiable par sa propriété d'absence de mémoire, ce qui correspond à beaucoup de situations concrètes. On s'intéresse au nombre $X_t$ de clients dans la file d'attente à l'instant $t$. Le processus $(X_t)_{t\ge0}$ est une chaîne de Markov à temps continu d'espace d'états $\mathbb{N}$. Dans une nomenclature due à Kendall, on dit qu'il s'agit d'une file M/M/1 de taux $\lambda$ et $\mu$: le premier M indique l'absence de mémoire<fnref n="2" /> des durées entre les arrivées, le second M indique l'absence de mémoire des durées de service, et le 1 final indique enfin qu'il n'y a qu'un seul serveur. Plus généralement, on peut définir des files d'attentes M/M/s, où s est un entier quelconque. Plus généralement encore la file M/M/s/K tient compte d'une taille maximale K de file d'attente ce qui modélise une salle d'attente de capacité limité. La file d'attente M/M/K/K modélise par exemple un parking avec K places. Il est possible de tenir compte de phénomènes supplémentaires comme l'impatience,

<footnote n="1">1. En anglais on dit « *queues* ».</footnote>

<footnote n="2">2. En anglais on dit « *Memoryless* ».</footnote>
la priorité, la politique de traitement 3, etc. De plus les files d'attentes peuvent être mises en tandem, voire en réseau, etc. Enfin, il est parfois plus réaliste de faire varier les intensités d'arrivée $\lambda$ et de service $\mu$ au cours du temps, ce qui mène mathématiquement à utiliser la notion naturelle de processus ponctuel.

Ce chapitre est consacré à la file d'attente M/M/∞, qui est sans doute la file la plus simple à étudier, et qui constitue un analogue à espace discret du processus d'Ornstein-Uhlenbeck. Dans le modèle de la file d'attente M/M/∞, les durées qui séparent les arrivées des clients sont i.i.d. de loi exponentielle de paramètre $\lambda$. Chaque client dispose d'un serveur dédié dès son arrivé qui débute son traitement immédiatement, et les durées de traitement sont i.i.d. de loi exponentielle de paramètre $\mu$, indépendantes du processus des arrivées. Une particularité de la file d'attente M/M/∞ est qu'à tout instant, le nombre de clients dans la file est égal au nombre de clients en cours de service.

La file d'attente M/M/∞ constitue une version simplifiée de la file d'attente M/M/K/K quand $\lambda/\mu$ est négligeable devant K. Elle modélise par exemple le nombre de voitures garées dans un parking de très grande capacité : les voitures sont les clients, les places du parking sont les serveurs, et la durée de service d'une voiture correspond à la durée d'occupation de la place de parking par la voiture. Des situations concrètes très variées peuvent bénéficier de la modélisation par la file d'attente M/M/∞.

Pour construire rigoureusement le processus $(X_t)_{t\ge0}$, on se place tout d'abord dans le cas où la file est vide au temps initial : $X_0 = 0$. Les durées séparant les arrivées des clients sont modélisées par une suite $(E_n)_{n\ge0}$ de v.a.r. i.i.d. de loi exponentielle de paramètre $\lambda$. Le $n$-client arrive donc au temps $T_n := E_1 + \dots + E_n$, qui suit une loi Gamma$(n, \lambda)$. Pour tout réel $t \ge 0$, le nombre de clients arrivés dans l'intervalle de temps $[0, t]$ est

$$N_t := \sum_{n=0}^{\infty} 1_{\{T_n \le t\}}.$$

Le processus de comptage $(N_t)_{t\ge0}$ est un processus de Poisson issu de 0 et d'intensité $\lambda$. Pour tout entier $n \ge 0$, on note $S_n$ la durée de service du $n$-client de sorte que ce client quitte la file au temps $E_n + S_n$. La suite $(S_n)_{n\ge0}$ est constituée de v.a.r. i.i.d. de loi exponentielle de paramètre $\mu$. On suppose que les variables aléatoires $X_0$, $(S_n)_{n\ge0}$, $(E_n)_{n\ge0}$ sont indépendantes. Pour tout $t \in \mathbb{R}_+$, le nombre de clients dans la file au temps $t$ vaut

$$\sum_{n=0}^{\infty} 1_{\{T_n \le t < T_n + S_n\}}.$$

Dans le cas général où $X_0$ n'est pas forcément nul, on introduit une suite $(S'_n)_{n\ge0}$ de v.a.r. i.i.d. de loi exponentielle de paramètre $\mu$, indépendante de $X_0$, et on modélise les durées de service de ces clients initiaux par $S'_0, \dots, S'_{X_0}$.

<footnote n="3">3. La politique de traitement standard est FIFO (pour «First In First Out»).</footnote>
On suppose que $X_0, (S_n)_{n \ge 0}, (S'_n)_{n \ge 0}, (E_n)_{n \ge 0}$ sont indépendantes. Pour tout réel $t \ge 0$, le nombre de clients dans la file au temps $t$ vaut donc

$$X_t := \sum_{n=0}^{X_0} 1_{\{S'_n > t\}} + \sum_{n=0}^{\infty} 1_{\{T_n \le t < T_n + S_n\}}.$$

## 11.1 Lois exponentielles

Les lemmes suivants sont très utiles, bien au-delà des files d'attentes. Ils jouent un rôle notamment dans la structure des processus de Markov.

**Lemme 11.1 (Absence de mémoire).** Si $E_1$ et $E_2$ sont deux variables aléatoires exponentielles indépendantes de paramètres respectifs $\lambda_1$ et $\lambda_2$, alors

$$\mathrm{Loi}(E_2 - E_1 | E_2 > E_1) = \mathrm{Exp}(\lambda_2) = \mathrm{Loi}(E_2)$$

c'est-à-dire que $\mathbb{P}(E_2 - E_1 > t | E_2 > E_1) = e^{-\lambda_2 t}$ pour tout $t \ge 0$.

*Démonstration.* Pour tout $t \ge 0$ on a $\mathbb{P}(E_2 > E_1 + t | E_2 > E_1) = e^{-\lambda_2 t}$ car

$$\mathbb{P}(E_2 > E_1) = \frac{\lambda_1}{\lambda_1 + \lambda_2}$$

tandis que $\mathbb{P}(E_2 > E_1 + t, E_2 > E_1) = \mathbb{P}(E_2 > E_1 + t)$ vaut

$$\iint_{v>u+t} \lambda_1 e^{-\lambda_1 u} \lambda_2 e^{-\lambda_2 v} dudv = \frac{\lambda_1}{\lambda_1 + \lambda_2} e^{-\lambda_2 t}.$$

□

**Lemme 11.2 (Horloges exponentielles en compétition).** Si $(E_i)_{i \in \mathcal{I}}$ est une famille finie ou infinie dénombrable de v.a.r. indépendantes de lois exponentielles de paramètres respectifs $(\lambda_i)_{i \in \mathcal{I}}$ avec $\lambda := \sum_{i \in \mathcal{I}} \lambda_i < \infty$, alors la variable aléatoire $M := \inf_{i \in \mathcal{I}} E_i$ suit la loi exponentielle $\mathrm{Exp}(\lambda)$, et de plus, presque sûrement, l'infimum est atteint en un unique entier aléatoire $I$, indépendant de $M$, dont la loi est donnée par $\mathbb{P}(I = i) = \lambda_i/\lambda$ pour tout $i \in \mathcal{I}$.

*Démonstration du lemme 11.2.* Pour tout $i \in \mathcal{I}$ et tout $x \in \mathbb{R}_+$,

$$\mathbb{P}(\min_{k \neq i}(E_k) > E_i \geq x) = \lambda_i \int_x^\infty \prod_{k \neq i} \mathbb{P}(E_k > y) e^{-\lambda_i y} dy = \frac{\lambda_i}{\lambda} e^{-\lambda x}.$$

C'est la formule pour $\mathbb{P}(I = i, M \ge x)$, qui donne la loi du couple $(I, M)$. □
**Lemme 11.3** (Lemme de Rényi sur la statistique d'ordre exponentielle). Si $E_1, \dots, E_n$ sont des v.a.r. i.i.d. de loi exponentielle Exp($\lambda$), et si

$$ \min(E_1, \dots, E_n) = E_{(1)} < \dots < E_{(n)} = \max(E_1, \dots, E_n) $$

est le réordonnement croissant (statistique d'ordre) de $E_1, \dots, E_n$, alors

$$ (E_{(1)}, \dots, E_{(n)}) \stackrel{\text{loi}}{=} \left( \frac{E_1}{n}, \dots, \frac{E_1}{n} + \dots + \frac{E_n}{n} \right). $$

En particulier, pour tout $1 \le k \le n$, avec la convention $E_0 := 0$,

$$ E_{(k)} \stackrel{\text{loi}}{=} \frac{E_1}{n} + \dots + \frac{E_k}{n-k+1} = E_{(k-1)} + \frac{E_k}{n-k+1}. $$

**Démonstration.** Pour le voir, on écrit, pour toute fonction mesurable et bornée $\varphi : \mathbb{R}^n \to \mathbb{R}$, en utilisant des changements de variables,

$$ \begin{align*} \mathrm{E}(\varphi(E_{(1)}, \dots, E_{(n)})) &= \int_{[0, \infty)^n} \varphi(x_{(1)}, \dots, x_{(n)}) \lambda^n e^{-\lambda(x_1+\dots+x_n)} dx_1 \dots dx_n \\ &= n! \int_{0 \le y_1 < \dots < y_n} \varphi(y_1, \dots, y_n) \lambda^n e^{-\lambda(y_1+\dots+y_n)} dy_1 \dots dy_n \\ &= \int_{[0, \infty)^n} \varphi(z_1, \dots, z_1 + \dots + z_n) n! \lambda^n e^{-\lambda(nz_1+\dots+1z_n)} dz, \end{align*} $$

où $dz := dz_1 \cdots dz_n$, ce qui signifie que les variables aléatoires

$$ E_{(1)}, E_{(2)} - E_{(1)}, \dots, E_{(n)} - E_{(n-1)} $$

sont indépendantes avec $E_{(k)} - E_{(k-1)} \sim \mathrm{Exp}((n-k+1)\lambda)$, $1 \le k \le n$. $\square$

Notons au passage que l'identité en loi concernant le minimum

$$ \min(E_1, \dots, E_n) = E_{(1)} \stackrel{\text{loi}}{=} \frac{1}{n} E_1 $$

est un cas spécial du lemme 11.2 avec $I = \{1, \dots, n\}$ et $\lambda_1 = \dots = \lambda_n$. D'autre part, l'identité en loi concernant le maximum

$$ \max(E_1, \dots, E_n) = E_{(n)} \stackrel{\text{loi}}{=} \frac{E_1}{n} + \dots + \frac{E_n}{1} $$

est démontrée autrement dans la preuve du théorème 10.8 (transformée de Laplace) ainsi que dans la preuve du théorème 13.4 (via la densité, par récurrence sur $n$). Elle apparaît également dans les chapitres 1, 9, et 17 (processus de Yule). Nous l'utilisons dans la section suivante pour le temps de demi-vie.

**Lemme 11.4** (Stabilité par sommation géométrique). Soient $G, E_1, E_2, \dots$ des variables aléatoires indépendantes avec $G$ de loi géométrique de paramètre $p \in ]0, 1]$ et $E_1, E_2, \dots$ de loi exponentielles de paramètre $\lambda$. Alors la somme aléatoire $E_1 + \dots + E_G$ suit la loi exponentielle de paramètre $p\lambda$.
Comme $E_1+\cdots+E_n \sim \text{Gamma}(n, \lambda)$, cela nous dit qu’un mélange géométrique de lois Gamma de même paramètre d’échelle est une loi exponentielle.

*Démonstration.* Pour tout $t \in \mathbb{R}$ et tout $k \ge 1$, on a $\mathbb{E}(e^{itE_k}) = \lambda/(\lambda - it)$, et

$$ 
\begin{align*}
\mathbb{E}(e^{it(E_1+\cdots+E_G)}) &= \sum_{n=1}^{\infty} \mathbb{E}(e^{it(E_1+\cdots+E_n)}) \mathbb{P}(G=n) \\
&= \sum_{n=1}^{\infty} \mathbb{E}(e^{itE_1})^n (1-p)^{n-1} p \\
&= \frac{p \mathbb{E}(e^{itE_1})}{1 - \mathbb{E}(e^{itE_1})(1-p)} \\
&= \frac{\lambda p}{\lambda p - it}.
\end{align*}
$$

□

Le lemme 11.4 indique que si on accepte chaque client du processus de Poisson des arrivées (de paramètre $\lambda$) avec probabilité $p$ alors le processus «aminci» résultant est un processus de Poisson d'intensité $p\lambda$.

Le lemme suivant achève notre petit catalogue de propriétés fondamentales des lois exponentielles. Il implique par exemple que les trajectoires du processus de Poisson $(N_t)_{t \ge 0}$ sont bien définies (elles n’explosent pas).

**Lemme 11.5 (Loi du zéro-un).** Si $(E_n)_{n \ge 1}$ sont des v.a.r. indépendantes de lois exponentielles de paramètres respectifs $(\lambda_n)_{n \ge 1}$ avec $0 < \lambda_n < \infty$ pour tout $n \ge 1$, alors l’événement $\{\sum_{n=1}^{\infty} E_n = \infty\}$ est de probabilité 0 ou 1 et

$$ \mathbb{P}\left(\sum_{n=1}^{\infty} E_n = \infty\right) = 1 \quad \text{si et seulement si} \quad \sum_{n=1}^{\infty} \frac{1}{\lambda_n} = \infty. $$

La série $\sum_{n=1}^{\infty} E_n$ diverge p.s. si et seulement si elle diverge en moyenne.

*Démonstration.* La variable $T_\infty = \sum_{n=1}^\infty E_n$ prend ses valeurs dans l'ensemble $\mathbb{R}_+ \cup \{\infty\}$. Le théorème de convergence monotone donne $\mathbb{E}(T_\infty) = \sum_{n=1}^\infty 1/\lambda_n$ dans $\mathbb{R}_+ \cup \{\infty\}$. Par conséquent, si $\sum_{n=1}^\infty 1/\lambda_n < \infty$ alors $\mathbb{E}(T_\infty) < \infty$ et donc $\mathbb{P}(T_\infty = \infty) = 0$. L'indépendance n'a pas été utilisée. Réciproquement, le théorème de convergence monotone et l'indépendance donnent

$$ \mathbb{E}[e^{-T_\infty}] = \prod_{n=1}^{\infty} \mathbb{E}[e^{-E_n}] = \prod_{n=1}^{\infty} \left(1 + \frac{1}{\lambda_n}\right)^{-1}. $$

Or $\sum_{n=1}^\infty 1/\lambda_n = \infty$ si et seulement si $\prod_{n=1}^\infty (1 + 1/\lambda_n)^{-1} = 0$, car pour toute suite $(a_n)_{n \ge 1}$ de réels positifs, le produit $\prod_{n=1}^\infty (1 + a_n)$ converge si et seulement si la série $\sum_{n=1}^\infty a_n$ converge. Ainsi, si $\sum_{n=1}^\infty 1/\lambda_n = \infty$, alors $\mathbb{E}(\exp(-T_\infty)) = 0$, et donc $\mathbb{P}(T_\infty = \infty) = 1$. □
## 11.2 Temps de demi-vie des clients initiaux

Si $\mu = 0$ alors le processus $(X_t)_{t \ge 0}$ coïncide avec le processus de Poisson $(X_0 + N_t)_{t \ge 0}$ issu de $X_0$ et d'intensité $\lambda$, où $(N_t)_{t \ge 0}$ est le processus de Poisson des arrivées, d'où $X_t \to +\infty$ p.s. quand $t \to \infty$ si $X_0 > 0$ et $\lambda > 0$.

On s'intéresse plutôt dans cette section au cas où $\lambda = 0$ et $X_0 = N$. La file comporte $N$ clients au temps $t = 0$, et aucun nouveau client n'arrive. Dans ce cas, le processus $(X_t)_{t \ge 0}$ est un processus de mort pur, et $X_t$ suit la loi $\text{Bin}(N, e^{-\mu t})$, d'où $X_t \to 0$ p.s. quand $t \to \infty$ si $\mu > 0$ ou si $X_0 = 0$. Avec les notations de l'introduction, les durées de service des $N$ clients initiaux sont $S'_1, \dots, S'_N$. Les instants de fin de service successifs sont notés $(T_n^{(N)})_{1 \le n \le N}$ :

$$ \min(S'_1, \dots, S'_n) = T_1^{(N)} \le \dots \le T_N^{(N)} = \max(S'_1, \dots, S'_n). $$

On a $X_t = 0$ si $t > T_N^{(N)}$. Le lemme de Rényi 11.3 indique que les v.a.r.

$$ (\tau_n^{(N)})_{1 \le n \le N} \quad \text{où} \quad \tau_n^{(N)} := T_n^{(N)} - T_{n-1}^{(N)} \quad \text{et} \quad T_0^{(N)} := 0 $$

sont indépendantes de lois exponentielles de paramètres respectifs

$$ ((N - n + 1)\mu)_{1 \le n \le N}. $$

On s'intéresse au temps de demi-vie $T_{\lfloor N/2 \rfloor}^{(N)}$, c'est-à-dire au temps au bout duquel la moitié des $N$ clients initiaux a quitté la file.

**Théorème 11.6** (Convergence du temps de demi-vie). On a

$$ T_{\lfloor N/2 \rfloor}^{(N)} \xrightarrow[N \to \infty]{\text{p.s.}} \frac{\log(2)}{\mu}. $$

*Démonstration.* Pour alléger les notations et simplifier, on omet l'exposant «(N)», on suppose que $N$ est pair, et on pose $M := N/2$. On observe tout d'abord que

$$ \begin{aligned} \mathbb{E}(T_{N/2}) &= \sum_{n=1}^{N/2} \mathbb{E}(\tau_n) = \frac{1}{\mu} \sum_{n=1}^{N/2} \frac{1}{N - n + 1} \\ &= \frac{1}{\mu} (H(N) - H(N/2)) = \frac{\log(2)}{\mu} + o(1), \end{aligned} $$

où $H(n)$ est la série harmonique. De plus,

$$ T_M - \mathbb{E}(T_M) = \sum_{n=1}^{M} (\tau_n - \mathbb{E}(\tau_n)), $$

où les variables $(\tau_n - \mathbb{E}(\tau_n))_{n \ge 1}$ sont indépendantes et centrées et $\tau_n - \mathbb{E}(\tau_n)$ a la même loi que $(E-1)/\mu(N+1-n)$ où $E$ suit la loi exponentielle de paramètre 1. Pour établir la convergence, nous allons utiliser l'inégalité de Markov
avec un moment d'ordre suffisamment élevé pour assurer la convergence. Le moment d'ordre 4 se contrôle comme suit :

$$
\begin{align*}
\mathrm{E}[(T_M - \mathrm{E}(T_M))^4] &= \sum_{n=1}^{M} \mathrm{E}[(\tau_n - \mathrm{E}(\tau_n))^4] \\
&\quad + 6 \sum_{1 \le m \ne n \le M} \mathrm{E}[(\tau_n - \mathrm{E}(\tau_n))^2] \mathrm{E}[(\tau_m - \mathrm{E}(\tau_m))^2] \\
&= \frac{\mathrm{E}[(E-1)^4]}{\mu^4} \sum_{n=M+1}^{N} \frac{1}{n^4} \\
&\quad + 6 \frac{(\mathrm{E}[(E-1)^2])^2}{\mu^4} \sum_{1+M \le m \ne n \le N} \frac{1}{n^2} \frac{1}{m^2} \\
&\le \frac{9}{\mu^4} \sum_{n=M+1}^{N} \frac{1}{n^4} + \frac{6}{\mu^4} \left( \sum_{n=M+1}^{N} \frac{1}{n^2} \right)^2,
\end{align*}
$$

puisque $E((E-1)^2) = 1$ et $E((E-1)^4) = 9$. Ainsi

$$
E[(T_M - \mathrm{E}(T_M))^4] = \mathcal{O}(N^{-2}).
$$

L'inégalité de Markov assure alors que, pour tout $r \ge 0$,

$$
P(|T_M - \mathrm{E}(T_M)| \ge r) \le \frac{C}{r^4 N^2}.
$$

Le lemme de Borel-Cantelli donne enfin $T_M \xrightarrow[N \to \infty]^{\text{p.s.}} \log(2)/\mu$. $\square$

**Remarque 11.7** (Demi-vie et radioactivé). La durée de vie des atomes de Carbone 14 est bien modélisée par une variable aléatoire de loi exponentielle de demi-vie 5734 ans. Ceci permet aux archéologues de dater les échantillons de matière organique.

**11.3 Loi du nombre de clients**

Cette section est dédiée au calcul de la loi de $X_t$. Ce calcul est mené avec deux méthodes différentes.

**Théorème 11.8** (Loi de $X_t$). Pour tout $k \in \mathbb{N}$ et tout $t \ge 0$,

$$
\mathrm{Loi}(X_t | X_0 = k) = \mathrm{Bin}(k, e^{-\mu t}) * \mathrm{Poi}\left(\frac{\lambda}{\mu}(1 - e^{-\mu t})\right)
$$

avec la convention $\mathrm{Bin}(0,p) = \delta_0$. De plus, si $X_0 \sim \mathrm{Poi}(\alpha)$ alors

$$
X_t \sim \mathrm{Poi}\left(\frac{\lambda}{\mu}(1 - e^{-\mu t}) + \alpha e^{-\mu t}\right).
$$
Démonstration. Adoptons les notations de l'introduction. Par indépendance,

$$
\text{Loi}(X_t | X_0 = k) = \text{Loi} \left( \sum_{n=0}^{X_0} 1_{\{S'_n > t\}} \mid X_0 = k \right) * \text{Loi} \left( \sum_{n=0}^{\infty} 1_{\{T_n \le t < T_{n+S_n}\}} \right).
$$

À nouveau par indépendance (de $X_0$ et de $(S'_n)_{n\ge 0}$ cette fois) on a

$$
\text{Loi}\left(\sum_{n=0}^{X_0} 1_{\{S'_n > t\}} \mid X_0 = k\right) = \text{Loi}\left(\sum_{n=0}^{k} 1_{\{S'_n > t\}}\right) = \text{Bin}(k, e^{-\mu t})
$$

car $1_{\{S'_1>t\}}, \dots, 1_{\{S'_k>t\}}$ sont des variables aléatoires i.i.d. de loi de Bernoulli
$\text{Ber}(e^{-\mu t})$. Il ne reste plus qu'à établir que

$$
\text{Loi}(Y_t) = \text{Poi}\left(\frac{\lambda}{\mu}(1 - e^{-\mu t})\right) \quad \text{où} \quad Y_t := \sum_{n=0}^{\infty} 1_{\{T_n \le t < T_{n+1}\}}.
$$

Rappelons que $(N_t)_{t \ge 0}$ est le processus de Poisson qui compte les temps $(T_n)_{n \ge 1}$. Pour $t > 0$, la variable aléatoire $N_t = \text{card}\{n \ge 1 : T_n \in [0, t]\}$ suit la loi de Poisson de paramètre $\lambda t$. De plus, sachant que $N_t = n$, la loi des $n$ instants de saut est celle d'un $n$-échantillon réordonné de v.a.r. i.i.d. de loi uniforme sur $[0, t]$. D'autre part, un client arrivé à un instant aléatoire de loi uniforme sur $[0, t]$ est encore présent dans la file à l'instant $t$ avec probabilité

$$
q(t) := \frac{1}{t} \int_{0}^{t} (1 - F_{\mu}(s)) \, ds = \frac{1}{\mu t} (1 - e^{-\mu t}),
$$

où $F_\mu$ est la fonction de répartition de la loi $\text{Exp}(\mu)$. On obtient donc, par
indépendance, pour tout $k \in \mathbb{N}$,

$$
\begin{align*}
\mathbb{P}(Y_t = k) &= \sum_{n=k}^{\infty} \mathbb{P}(Y_t = k | N_t = n) \mathbb{P}(N_t = n) \\
&= \sum_{n=k}^{\infty} \binom{n}{k} q(t)^k (1-q(t))^{n-k} e^{-\lambda t} \frac{(\lambda t)^n}{n!} \\
&= e^{-\lambda t q(t)} \frac{(\lambda t q(t))^k}{k!},
\end{align*}
$$

ce qui montre bien que $Y_t$ suit la loi de Poisson de paramètre $\lambda(1-e^{-\mu t})/\mu$. $\square$

**Remarque 11.9** (Loi de service quelconque et file d'attente M/G/∞). Si les durées de service sont i.i.d. de loi de fonction de répartition F quelconque, alors la preuve du théorème (11.8) permet d'établir que la loi de $X_t$ sachant $X_0 = 0$ est la loi de Poisson de paramètre

$$
\lambda \int_{0}^{t} (1 - F(s)) \, ds.
$$

On parle dans ce cas de file d'attente M/G/∞, où G signifie «General».
Il s'avère que le théorème 11.8 peut également être obtenu en utilisant les équations de Chapman-Kolmogorov. Cette approche reste valable pour les chaînes de Markov à temps continu, bien au-delà du cas particulier étudié.

**Théorème 11.10** (Équations de Chapman-Kolmogorov). Pour tous $t \in \mathbb{R}_+$ et $n \in \mathbb{N}$, soit $p_n(t) := P(X_t = n)$. Les fonctions $(p_n)_{n \in \mathbb{N}}$ sont de classe $C^1$ et vérifient le système infini d'équations différentielles linéaires suivant :

$$ \forall n \in \mathbb{N}, \quad p'_{n}(t) = -(\lambda + n\mu)p_{n}(t) + \lambda p_{n-1}(t) + (n+1)\mu p_{n+1}(t), $$

avec $p_{-1}(t) := 0$ pour tout $t \ge 0$, la loi initiale $(p_n(0))_{n \in \mathbb{N}}$ étant donnée.

**Démonstration.** Pour calculer $p_n(t+h)$, on remarque que si $X_{t+h} = n$ alors l'une des conditions incompatibles suivantes est réalisée :

1. $X_u = n$ pour tout $u \in [t, t+h]$;

2. $X_t = n-1$ et une seule transition a lieu ($n-1 \to n$) dans l'intervalle de temps $[t, t+h]$;

3. $X_t = n+1$ et une seule transition a lieu ($n+1 \to n$) dans l'intervalle de temps $[t, t+h]$;

4. dans l'intervalle de temps $[t, t+h]$, au moins deux transitions ont lieu.

D'après les lemmes 11.1 et 11.2, le processus reste dans l'état $n$ un temps aléatoire de loi exponentielle de paramètre $\lambda + n\mu$ puis saute de $n$ à $n-1$ ou $n+1$ avec probabilités respectives $n\mu/(\lambda+n\mu)$ et $\lambda/(\lambda+n\mu)$. En conditionnant par l'évènement $\{X_t = n\}$, on en déduit l'égalité

$$ p_n(t+h) = p_n(t)\{1 - \lambda h - n\mu h\} + \lambda h p_{n-1}(t) + (n+1)\mu h p_{n+1}(t) + o(h). $$

Le résultat désiré s'obtient en faisant tendre $h$ vers 0. □

Le corollaire ci-dessous du théorème 11.10 affirme que la fonction généra-trice de $X_t$ est solution d'une équation aux dérivées partielles. Il fournit donc une nouvelle preuve du théorème 11.8, analytique, basée sur les fonctions géné-ratrices (ou transformée de Laplace). Pour tous $t \in \mathbb{R}_+$ et $s \in [-1, 1]$, soit

$$ G(s,t) := \mathbb{E}(s^{X_t}) = \sum_{n=0}^{\infty} p_n(t)s^n. $$

**Corollaire 11.11** (Fonction génératrice de $X_t$). Si $G_0$ est la fonction géné-ratrice de $X_0$ alors, pour tout $t \ge 0$ et tout $s \in [-1, 1]$,

$$ G(s,t) = G_0(1 - e^{-\mu t} + s e^{-\mu t}) \exp \left( \frac{\lambda}{\mu} (1 - e^{-\mu t}) (s - 1) \right). $$

**Démonstration.** Par les équations de Chapman-Kolmogorov (théorème 11.10),

$$ \partial_t G(s,t) = \sum_{n=0}^{\infty} p'_{n}(t)s^n $$
$$
\begin{align*}
&= \sum_{n=0}^{\infty} - (\lambda + n\mu)p_n(t) + \lambda p_{n-1}(t) + (n+1)\mu p_{n+1}(t))s^n \\
&= -\lambda(1-s) \sum_{n=0}^{\infty} p_n(t)s^n + \mu(1-s) \sum_{n=0}^{\infty} np_n(t)s^{n-1}.
\end{align*}
$$

La fonction G est donc solution de l'équation aux dérivées partielles suivante :

$$
\partial_t G(s,t) = (1-s)[\mu \partial_s G(s,t) - \lambda G(s,t)].
$$

On pose $H = \log(G)$ puis on effectue ensuite le changement de variables
$(s, \tau) = (s, (1-s)e^{-\mu t})$. En notant $J$ la fonction telle que $H(s,t) = J(s, \tau)$,
la fonction $J$ s'écrit sous la forme

$$
J(s, \tau) = \frac{\lambda}{\mu} s + h(\tau),
$$

où $h$ reste à déterminer. Ceci revient à dire que $H$ est de la forme

$$
H(s,t) = \frac{\lambda}{\mu}s + h((1-s)e^{-\mu t}).
$$

On conclut en remarquant que la fonction $h$ est caractérisée par le fait que
$H(s, 0)$ est égal à $\log(G_0(s))$ pour tout $s \in [-1, 1]$. $\square$

Le processus $(X_t)_{t\ge0}$ est une chaîne de Markov à temps continu. Le semi-
groupe markovien associé est la famille d'opérateurs $(P_t)_{t\ge0}$ définie comme
suit : pour tous $t \ge 0$ et $f : \mathbb{N} \to \mathbb{R}$, $P_t f$ est la fonction $\mathbb{N} \to \mathbb{R}$ donnée par

$$
P_t f(n) := \mathbb{E}(f(X_t) | X_0 = n).
$$

Si $\nu$ est une mesure de probabilité sur $\mathbb{N}$ on notera $\nu P_t$ la loi de $X_t$ sachant que $X_0$ suit la loi $\nu$. On a, pour toute fonction $f$ bornée,

$$
\nu P_t f = \sum_{n \in \mathbb{N}} \nu(n) P_t f(n).
$$

D'après le théorème 11.8, la mesure de probabilité $P_t(\cdot)(n)$ est la loi

$$
\text{Bin}(n, e^{-\mu t}) * \text{Poi}\left(\frac{\lambda}{\mu}(1 - e^{-\mu t})\right).
$$

**Remarque 11.12** (Propriétés fondamentales du semi-groupe markovien). *Le semi-groupe* $(P_t)_{t \ge 0}$ *vérifie les propriétés fondamentales suivantes :*

1. Élément neutre : $P_0 = I$ c'est-à-dire que $P_0 f = f$ pour tout $f$ ;

2. Positivité et normalisation : si $f \ge 0$ alors $P_t f \ge 0$, et $P_t 1 = 1$;

3. Propriété de semi-groupe : pour toute fonction $f$ et tous $s, t \ge 0$,

$P_{t+s}f = (P_t \circ P_s)f;$
4. Équation de Chapman-Kolmogorov : pour toute fonction $f$,

$$ \partial_t P_t f = LP_t f = P_t L f; $$

où $L$ est appelé générateur infinitésimal de $X$ et est défini par

$$ Lf(n) := \lambda(f(n+1) - f(n)) + n\mu(f(n-1) - f(n)). $$

En fait, formellement, tout se passe comme si $P_t = e^{tL}$. Les deux premiers points sont évidents. Le troisième découle de la formule explicite de $P_t$. C'est une reformulation de la propriété de Markov. La dernière assertion s'obtient grâce à 3. et un raisonnement analogue à celui de la preuve du théorème 11.10.

**Remarque 11.13** (Lire la dynamique dans le générateur infinitésimal). Si on définit le noyau de transition

$$ Q(n, \cdot) := \frac{\lambda}{\lambda + n\mu} \delta_{n+1} + \frac{n\mu}{\lambda + n\mu} \delta_{n-1} $$

et la «temporisation» $D(n) := \lambda + n\mu$ alors le générateur infinitésimal s'écrit

$$ Lf(n) = D(n) \int (f(m) - f(n)) Q(n, dm). $$

On lit dans le générateur infinitésimal un algorithme de simulation des trajec-
toires du processus : sachant qu'il est en $n$, le processus $X$ saute au bout d'un
temps exponentiel de paramètre $D(n)$ vers $n+1$ avec probabilité $\lambda/(\lambda + n\mu)$
et vers $n-1$ avec probabilité $n\mu/(\lambda + n\mu)$. On trouvera par exemple dans le
chapitre 9 page 120 un autre exemple de ce type. Notons qu'on peut voir $L$
comme une matrice avec une infinité de lignes et une infinité de colonnes, en
posant $L(n, m) := L1_n(m)$ de sorte que $Lf(n) = \sum_{m \in \mathbb{N}} L(n, m)f(m)$.

On peut associer à tout processus de Markov raisonnable son semi-groupe
et son générateur infinitésimal. Si ce dernier est souvent explicite, il est rare
que ce soit le cas pour les mesures $P_t(\cdot)(x)$. En ce sens, le processus étudié
dans ce chapitre est assez remarquable.

## 11.4 Comportement en temps long

On adopte les notations suivantes :

$$ \rho := \frac{\lambda}{\mu} \quad \text{et} \quad \pi := \text{Poi}(\rho). $$

Si $X_0 \sim \pi$, alors le théorème 11.8 assure que $X_t \sim \pi$ pour tout $t \ge 0$. Ainsi
la loi $\pi$ est invariante sous l'action de $P_t$. De manière équivalente, $\pi L = 0$.
Comme pour tous les processus de naissance et mort, la loi invariante $\pi$, si elle existe, est de plus réversible pour $L$ (ou $P_t$), c'est-à-dire que $L$ et $P_t$ sont auto-adjoints dans $L^2(\pi)$. Enfin, pour toute loi initiale (de $X_0$),

$$X_t \xrightarrow{\text{loi}} \pi.$$

On souhaite ici quantifier cette convergence en loi.

Soit $\mathcal{P}_1(\mathbb{N})$ l'ensemble des mesures de probabilité $\nu$ sur $\mathbb{N}$ possédant un moment d'ordre 1, c'est-à-dire l'ensemble des $\nu : \mathbb{N} \to [0, 1]$ vérifiant $\sum_{x \in \mathbb{N}} \nu(x) = 1$ et $\sum_{x \in \mathbb{N}} x\nu(x) < \infty$. La distance de Wasserstein sur $\mathcal{P}_1(\mathbb{N})$ est définie pour tous $\nu, \tilde{\nu} \in \mathcal{P}_1(\mathbb{N})$ par la formule variationnelle de couplage

$$W_1(\nu, \tilde{\nu}) := \inf \{ \mathbb{E}|X - \tilde{X}| : X \sim \nu, \tilde{X} \sim \tilde{\nu} \}.$$

Relions cette distance à la distance en variation totale introduite dans le chapitre 1. Pour des v.a. entières, $X \neq \tilde{X}$ implique $|X - \tilde{X}| \ge 1$. On a donc

$$d_{\text{VT}}(\nu, \tilde{\nu}) \le W_1(\nu, \tilde{\nu}).$$

**Théorème 11.14** (Contraction de Wasserstein). Si $\nu, \tilde{\nu} \in \mathcal{P}_1(\mathbb{N})$ alors

$$W_1(\nu P_t, \tilde{\nu} P_t) \le e^{-\mu t} W_1(\nu, \tilde{\nu}),$$

où, pour tout $\eta \in \mathcal{P}(\mathbb{N})$, $\eta P_t$ est la loi de $X_t$ si $X_0$ suit la loi $\eta$.

**Démonstration.** On procède par couplage. Supposons dans un premier temps que $\nu = \delta_m$ et $\tilde{\nu} = \delta_n$ avec $m < n$. Soit $X_t$ de loi $\delta_m P_t$ et $B_t$ de loi $\text{Bin}(n-m, e^{-\mu t})$ indépendante de $X_t$. D'après le théorème 11.8, la variable aléatoire $\tilde{X}_t = X_t + B_t$ suit la loi $\delta_n P_t$. On a donc

$$\mathbb{E}|X_t - \tilde{X}_t| = \mathbb{E}(B_t) = (n-m)e^{-\mu t}.$$

Comme $(X_t, \tilde{X}_t)$ est un couplage de $(\delta_m P_t, \delta_n P_t)$, on a

$$W_1(\delta_m P_t, \delta_n P_t) \le e^{-\mu t} (n-m) = e^{-\mu t} W_1(\delta_m, \delta_n).$$

Concluons à présent dans le cas général. Soit $\nu$ et $\tilde{\nu}$ dans $\mathcal{P}_1(\mathbb{N})$. Soit $(X_0, \tilde{X}_0)$ un couplage de $\nu$ et $\tilde{\nu}$. Comme ci-dessus, on construit $(X_t)$ et $(\tilde{X}_t)$ tels que

$$\text{Loi}(|X_t - \tilde{X}_t| | X_0, \tilde{X}_0) = \text{Bin}(|X_0 - \tilde{X}_0|, e^{-\mu t}).$$

On en déduit donc

$$W_1(\nu P_t, \tilde{\nu} P_t) \le e^{-\mu t} \mathbb{E}|X_0 - \tilde{X}_0|.$$

L'infimum sur tous les couplages de $X_0$ et $\tilde{X}_0$ donne la propriété attendue. $\square$

Ce résultat de contraction fournit une estimation pour la convergence, au sens de la distance de Wasserstein, de la loi de $X_t$ vers la loi invariante $\pi$ de $X$. La vitesse de convergence obtenue ne dépend pas du taux de naissance $\lambda!$
**Corollaire 11.15 (Convergence à l'équilibre).** Si $\nu \in \mathcal{P}_1(\mathbb{N})$ alors

$$W_1(\nu P_t, \pi) \le e^{-\mu t} W_1(\nu, \pi).$$

**Remarque 11.16 (Optimalité de la vitesse de convergence).** La constante $\mu$ dans la convergence à l'équilibre à vitesse exponentielle $e^{-\mu t}$ est optimale. En effet, supposons que la convergence à l'équilibre ait lieu avec une constante $\mu' > \mu$. Alors, en prenant pour loi initiale $\nu = \delta_n$ avec $n$ supérieur à la moyenne $\rho$ de la loi invariante $\pi$, et en notant $X$ une v.a. de loi $\pi$, on aurait

$$(n - \rho)e^{-\mu t} = \mathbb{E}_n(X_t) - \rho \le W_1(\delta_n P_t, \pi) \le e^{-\mu' t}W_1(\nu, \pi) = \mathbb{E}|n - X|e^{-\mu' t},$$

ce qui conduit à une contradiction car $n - \rho > 0$ et $e^{-\mu' t}/e^{-\mu t} \xrightarrow{t \to \infty} 0$.

Dans la même veine, on peut montrer que la fonction $P_t f$ converge vers la fonction constante égale à l'intégrale de $f$ pour $\pi$ en montrant que son gradient (discret) converge vers 0.

**Théorème 11.17 (Commutation).** Si $f : \mathbb{N} \to \mathbb{R}$ alors pour tout $t \ge 0$,

$$DP_t f(n) = e^{-\mu t} P_t (Df)(n),$$

où $Dg$ est la fonction définie sur $\mathbb{N}$ par $Dg(n) = g(n + 1) - g(n)$.

**Démonstration.** D'après le théorème 11.8,

$$P_t f(n + 1) = \mathbb{E}f(Z_1 + \cdots + Z_{n+1} + Y),$$

où $Z_1, \dots, Z_{n+1}$ sont des variables aléatoires indépendantes de loi $Ber(e^{-\mu t})$ et indépendantes de $Y$ de loi de Poisson de paramètre $\rho(1 - e^{-\mu t})$. On a donc

$$
\begin{align*}
P_t f(n+1) &= \mathbb{E}(f(Z_1 + \cdots + Z_{n+1} + Y)) \\
&= e^{-\mu t} \mathbb{E}(f(Z_1 + \cdots + Z_n + 1 + Y)) \\
&\quad + (1 - e^{-\mu t}) \mathbb{E}(f(Z_1 + \cdots + Z_n + Y)) \\
&= e^{-\mu t} P_t(Df)(n) + P_t f(n),
\end{align*}
$$

ce qui fournit le résultat. $\square$

Si $f$ est une fonction 1-lipschitzienne, au sens où $|Df(n)| \le 1$ pour tout $n \in \mathbb{N}$ alors le théorème précédent assure que $P_t f$ est une fonction $e^{-\mu t}$-lipschitzienne. On retrouve ainsi la convergence pour la distance de Wasserstein par la formulation duale de cette distance, de Kantorovitch-Rubinstein

$$W_1(\nu, \tilde{\nu}) = \sup \left\{ \iint f d\nu - \iint f d\tilde{\nu} : f \text{ 1-lipschitzienne} \right\}.$$
**Remarque 11.18** (Files d'attente, naissance et mort). La figure 11.1 compare les trajectoires des files M/M/1 et M/M/∞. Ces files d'attente font partie des processus de vie et de mort<fnref n="4" />, qui sont les processus de Markov à temps continu et d'espace d'états $\mathbb{N}$ à transitions aux plus proches voisins. Leur générateur infinitésimal a la forme suivante (avec la convention $\mu_0 = 0$) :

$$Lf(n) = \lambda_n(f(n+1) - f(n)) + \mu_n(f(n-1) - f(n)).$$

Le cas M/M/1 correspond à $n \in \mathbb{N} \mapsto \lambda_n$ et $n \in \mathbb{N}^* \mapsto \mu_n$ constantes, tandis que le cas M/M/∞ correspond à $n \in \mathbb{N} \mapsto \lambda_n$ constante et $\mu_n = n\mu$ pour tout $n \in \mathbb{N}$. Le cas où $n \mapsto \lambda_n$ est linéaire et $\mu \equiv 0$ correspond au processus de Yule (processus de Galton-Watson à temps continu de loi de reproduction $\delta_2$).

<figure><img src="image_35.png" /><figcaption>**Fig. 11.1.** Un début de trajectoire de file d'attente M/M/1 de paramètres 1 et 3, et de file d'attente M/M/∞ de paramètres 1 et .5. Les deux trajectoires sont issues de 100. Le premier processus a pour loi invariante la loi géométrique sur $\mathbb{N}$ de paramètre $1/3$, et le second la loi de Poisson de paramètre $1/.5 = 2$. Dans les deux cas, la loi invariante a pour moyenne 2. La première trajectoire à un comportement linéaire et la seconde un comportement exponentiel.</figcaption></figure>

<footnote n="4">4. On dit aussi «de naissance et de mort», «Birth and death» en anglais.</footnote>
## 11.5 Pour aller plus loin

Les résultats classiques concernant les lois exponentielles sont disponibles par exemple dans le livre de James Norris [Nor98a]. Dans le livre de Geoffrey Grimmett et David Stirzaker [GS01], une version à temps discret du processus X est étudiée. Pour l'anecdote, elle peut modéliser par exemple le nombre de coquilles présentes dans un manuscrit après *n* relectures : à chaque lecture, l'auteur corrige chaque coquille avec probabilité *p*, indépendamment des autres, et en ajoute par mégarde un nombre aléatoire de loi de Poisson (loi des événements rares ou loi des petits nombres).

Convenablement renormalisé, le processus à espace d'états discret converge vers le processus d'Ornstein-Uhlenbeck étudié dans le chapitre 25. On renvoie à ce sujet par exemple au livre de Philippe Robert [Rob00]. Ce livre propose plus généralement une étude accessible et assez générale des files d'attentes, au moyen de la théorie des martingales et des processus ponctuels. On peut également consulter le livre de François Baccelli et Pierre Brémaud [BB03].

On peut établir de nombreux autres résultats pour le processus X. On trouve par exemple dans [Cha06] des inégalités fonctionnelles liées au comportement en temps long. La relation de commutation du théorème 11.17 est une propriété remarquable du semi-groupe. Le processus d'Ornstein-Uhlenbeck satisfait la même relation en remplaçant le gradient discret D par un gradient usuel sur $\mathbb{R}$, et cette liaison est une instance d'un phénomène plus général qui fait l'objet du chapitre 27. On peut obtenir des versions un peu plus faibles pour des processus plus généraux qui fournissent des estimations pour la convergence vers la loi invariante, voir par exemple l'article de Pietro Caputo, Paolo Dai Pra, et Gustavo Posta [CDPP09].
<figure><img src="image_36.png" /></figure>
# Modèle de Wright-Fisher

**Mots-clés.** Dynamique de population ; génétique des population.

**Outils.** Chaîne de Markov ; martingale ; théorème d'arrêt ; décomposition spectrale.

**Difficulté.** **

Ce chapitre est consacré à des modèles d'évolution de génotypes au fil des générations dans une population de taille constante. Le point de départ est la loi de Hardy-Weinberg liée à la théorie de Mendel en population infinie. Nous présentons ensuite quelques modèles classiques en population finie : le modèle de Moran puis les modèles de Wright-Fisher. La dernière section est consacrée au modèle de Cannings qui est une généralisation des précédents.

## 12.1 Loi de Hardy-Weinberg

La théorie de Mendel assure la préservation de la variabilité des génotypes.

**Théorème 12.1** (Loi de Hardy-Weinberg pour la théorie de Mendel). Considérons les fréquences alléliques dans une population pour un gène pouvant s'exprimer sous la forme de deux allèles A et B dans une population hermaphrodite<fnref n="1" /> diploïde<fnref n="2" /> idéale :

1. les générations sont séparées (elles ne coexistent pas);

2. la population est de taille infinie;

3. les individus s'y unissent aléatoirement (pas de choix du conjoint);

4. il n'y a pas de migration, pas de mutation, pas de sélection.

<footnote n="1">1. Hermaphrodite : chaque individu peut être mâle ou femelle (bisexualité).</footnote>

<footnote n="2">2. Diploïde : qui possède un double assortiment de chromosomes semblables.</footnote>
Si $p \in [0, 1]$ est la proportion d'allèles A dans la population initiale et $1-p$ la proportion d'allèles B alors c'est encore le cas à chaque génération suivante. De plus, si $r, s, t$ sont les proportions respectives des génotypes AA, AB et BB dans la population initiale avec $r+s+t=1$ alors les fréquences de AA, AB, et BB pour toutes les générations suivantes sont égales à

$$ (r + \frac{s}{2})^2, \quad 2(r + \frac{s}{2})(t + \frac{s}{2}), \quad (t + \frac{s}{2})^2. $$

Notons que $p = r + s/2$ (contributions de AA et AB normalisées à 1).

*Démonstration.* Un individu AA est issu

— à coup sûr d'un couple AA et AA qui a une probabilité $r^2$ de se former,

— avec probabilité $1/2$ d'un couple AA et AB qui a une probabilité $2rs$ de se former,

— avec probabilité $1/4$ d'un couple AB et AB qui a une probabilité $s^2$ de se former.

La fréquence du génotype AA à la génération 1 est donc donnée par

$$ r^2 + rs + \frac{s^2}{4} = \left(r + \frac{s}{2}\right)^2. $$

Par symétrie, la fréquence du génotype BB est $(t + s/2)^2$. La dernière est obtenue par différence. En réappliquant ces formules on voit que les proportions sont constantes dès la première génération. La fréquence de l'allèle A dans la population initiale est $2r + s$ (normalisée à 2 car chaque individu a deux allèles). À la génération suivante, elle vaut

$$ 2\left(r + \frac{s}{2}\right)^2 + 2\left(r + \frac{s}{2}\right)\left(t + \frac{s}{2}\right) = 2r + s. $$

Elle est donc bien constante au cours du temps. □

La loi de Hardy-Weinberg est généralisable à un locus<fnref n="3" /> avec plusieurs allèles $A_1, \dots, A_k$. Les relations de dominance entre allèles n'ont pas d'effet sur l'évolution des fréquences alléliques. Dans la pratique, on observe le phénoène de Hardy-Weinberg dans des populations variées au-delà des hypothèses que nous avons faites ici. Il faut beaucoup s'éloigner de ces hypothèses pour briser cet équilibre et obtenir une évolution au cours des générations. C'est ce qui motive l'introduction des modèles étudiés dans la suite du chapitre, qui décrivent l'évolution d'une population de taille finie.

## 12.2 Modèle de Moran

Le modèle de Moran est un modèle extrêmement simple de l'évolution de la fréquence de différents allèles d'un même gène au cours du temps. On

<footnote n="3">3. Locus : localisation précise d'un gène particulier sur un chromosome.</footnote>
considère dans cette section le cas de deux allèles A et B et on suppose que la population est de taille finie constante N, hermaphrodite, que les générations sont séparées, qu'il n'y a pas de migration, mutation, sélection. La transition de la génération n à la génération n + 1 se fait en tirant avec remise deux individus dans la génération n, le premier a deux fils qui reçoivent tous les deux l'allèle de leur père, le second n'a pas de descendance, et tous les autres individus ont exactement un seul descendant. On note $X_n$ le nombre d'allèle A dans la population à la génération n. Le nombre d'allèles A

— augmente d'une unité si le tirage avec remise donne (A, B),

— diminue d'une unité si le tirage est (B, A).

— reste constant sinon.

On peut réaliser $(X_n)_{n \ge 0}$ comme une suite récurrente aléatoire de la ma-
nière suivante. Soit $(U_n)_{n \ge 1}$ et $(V_n)_{n \ge 1}$ deux suites indépendantes de v.a.r.
i.i.d. de loi uniforme sur $[0, 1]$. Soit $X_0$ une v.a. à valeurs dans $\{0, 1, \dots, N\}$
indépendante des v.a. introduites ci-dessus. On pose

$$
X_{n+1} = X_n + \mathbb{1}_{\{U_{n+1} < X_n/N\}} - \mathbb{1}_{\{V_{n+1} < X_n/N\}}.
$$

La suite $(X_n)_{n \ge 0}$ est une chaîne de Markov sur $\{0, 1, \dots, N\}$ de matrice de transition $\mathbf{P}$ donnée pour tous $x, y \in \{0, 1, \dots, N\}$ par

$$
\mathbf{P}(x, y) = \frac{1}{N^2} \begin{cases} x(N-x) & \text{si } |x-y| = 1; \\ x^2 + (N-x)^2 & \text{si } y = x; \\ 0 & \text{sinon.} \end{cases}
$$

Les états 0 et N sont absorbants. Les autres mènent à {0, N} et sont donc
transitoires. Le premier instant où un seul allèle est présent est noté

$$
T := \inf\{n \ge 0 : X_n \in \{0, N\}\}.
$$

On dit que $T$ est le temps de fixation. L'espace d'états étant fini, le temps d'atteinte $T$ de $\{0, N\}$ est presque sûrement fini et intégrable :

$$
\mathbb{P}(T < \infty) = 1 \quad \text{et} \quad \mathbb{E}(T) < \infty.
$$

De plus $X_T$ suit une loi de Bernoulli portée par 0 et N et

$$
X_n \xrightarrow{n \to \infty} X_T.
$$

On a $\mathbb{P}(X_T = 0) + \mathbb{P}(X_T = N) = 1$. L'événement $\{X_T = N\}$ (respectivement $\{X_T = 0\}$) signifie que l'allèle A (respectivement B) est fixé.

Adoptons la notation standard $\mathbb{P}_x = \mathbb{P}(\cdot | X_0 = x)$ et $\mathbb{E}_x = \mathbb{E}(\cdot | X_0 = x)$.

**Théorème 12.2** (Lieu et temps de fixation). *Pour tout* $x \in \{0, 1, \dots, N\}$,

$$
\mathbb{P}_x(X_T = N) = 1 - \mathbb{P}_x(X_T = 0) = \frac{x}{N}.
$$

De plus, pour tout $p \in ]0, 1[$,

$$
\mathbb{E}_{\lfloor pN \rfloor}(T) \underset{N \to \infty}{\sim} -N^2(p \log p + (1-p) \log(1-p)).
$$
*Démonstration.* Plaçons-nous dans le cas où $X_0 \sim \delta_x$. On a $\mathbb{P}_x(T < \infty) = 1$, et comme $T$ est un temps de fixation, on sait que $(X_n)_{n \ge 0}$ converge p.s. vers $X_T$. Comme l'espace d'états est borné, le théorème de convergence dominée s'applique et la convergence a lieu dans $L^1$, d'où

$$ \lim_{n \to \infty} \mathbb{E}_x(X_n) = \mathbb{E}(X_T) = 0 \mathbb{P}_x(X_T = 0) + N \mathbb{P}_x(X_T = N) = N \mathbb{P}_x(X_T = N). $$

Or la formule pour $\mathbf{P}$ donne, pour tout $n \ge 1$,

$$ \mathbb{E}(X_n | X_{n-1} = x) = \frac{(x(N-x)(x-1+x+1) + (x^2 + (N-x)^2)x)}{N^2} = x, $$

d'où

$$ \mathbb{E}(X_n) = \mathbb{E}(\mathbb{E}(X_n | X_{n-1})) = \mathbb{E}(X_{n-1}) = \dots = \mathbb{E}(X_0) = x, $$

ce qui donne enfin $\mathbb{P}_x(X_T = N) = x/N = 1 - \mathbb{P}_x(X_T = 0)$. Notons que ces formules peuvent aussi s'obtenir en utilisant le système linéaire vérifié par le vecteur des probabilités d'atteinte pour la chaîne de Markov $(X_n)_{n \ge 0}$.

Notons $m(x) = \mathbb{E}_x(T)$ le temps moyen de fixation pour la chaîne issue de $x$. On a $m(0) = m(N) = 0$ et $m(x) = 1 + \sum_{y=0}^N \mathbf{P}(x, y)m(y)$ pour tout $x \in \{1, \dots, N-1\}$. On obtient donc une formule de récurrence à trois termes :

$$ m(x+1) - 2m(x) + m(x-1) = - \frac{N^2}{x(N-x)}. $$

On en déduit que

$$ m(x) = N \left( \sum_{y=1}^{x} \frac{N-x}{N-y} + \sum_{y=x+1}^{N-1} \frac{x}{y} \right). $$

Soit $0 < p < 1$. Pour tout $N \ge 1$, posons $x_N = \lfloor pN \rfloor$. On a alors

$$ \begin{align*} \frac{m(x_N)}{N^2} &= \frac{1}{N} \sum_{y=1}^{x_N} \frac{N-x_N}{N-y} + \frac{1}{N} \sum_{y=x_N+1}^{N-1} \frac{x_N}{y} \\ &= \frac{1 - N^{-1}x_N}{N} \sum_{y=1}^{x_N} \frac{1}{1 - N^{-1}y} + \frac{N^{-1}x_N}{N} \sum_{y=x_N+1}^{N-1} \frac{1}{N^{-1}y} \\ &\underset{N \to \infty}{\longrightarrow} (1-p) \underbrace{\int_0^p \frac{1}{1-t} dt}_{=-(1-p) \log(1-p) + p \log(p)} + p \int_p^1 \frac{1}{t} dt. \end{align*} $$

Ainsi, dans une population de taille $N$ avec $N$ grand, le temps moyen d'absorption partant de $x = pN$ avec $0 < p < 1$ est de l'ordre de

$$ -N^2(p \log p + (1-p) \log(1-p)). $$

□
La fonction $p \in [0, 1] \mapsto -(p \log(p) + (1-p) \log(1-p))$ est l'entropie de Boltzmann ou de Shannon de la loi de Bernoulli $\text{Ber}(p)$. Elle est continue, strictement concave, positive, symétrique par rapport à $1/2$, et atteint son maximum $\log(2)$ pour $p = 1/2$, et son minimum 0 pour $p = 0$ et $p = 1$.

**Remarque 12.3 (Martingale).** Soit $\mathcal{F}_n = \sigma(X_0, \dots, X_n)$ pour tout $n \ge 0$. En reprenant la preuve du théorème 12.2, on a, pour tout $n \ge 1$,

$$ \mathbb{E}(X_n | \mathcal{F}_{n-1}) = \mathbb{E}(X_n | X_{n-1}) = X_{n-1}. $$

Ainsi, en plus d'être une chaîne de Markov, la suite $(X_n)_{n \ge 0}$ est une martingale pour sa filtration naturelle $(\mathcal{F}_n)_{n \ge 0}$. En particulier, on a bien sûr la conservation $\mathbb{E}(X_0) = \mathbb{E}(X_n)$ pour tout $n \ge 1$. Le théorème d'arrêt indique que $(X_{n \wedge T})_{n \ge 0}$ est aussi une martingale et donc

$$ \mathbb{E}(X_0) = \mathbb{E}(X_{n \wedge T}), $$

qui converge vers $\mathbb{E}(X_T)$ par convergence dominée quand $n \to \infty$ car $T$ est fini p.s. et $X$ est uniformément bornée. Comme $T$ est un temps de fixation, on a en fait $X_n = X_{n \wedge T}$. Il est également possible d'appliquer le théorème de convergence p.s. et dans $L^p$ des martingales bornées dans $L^p$ avec $p > 1$.

## 12.3 Modèle de Wright-Fisher et fixation

Le modèle de Wright-Fisher est un modèle simple de l'évolution de la fréquence d'un gène à plusieurs allèles au fil des générations dans une population. On suppose que la population est de taille finie constante, hermaphrodite, que les générations sont séparées, qu'il n'y a pas de migration et que les unions sont indépendantes du caractère étudié.

On néglige dans un premier temps les phénomènes de mutation et sélection. On note $N \in \mathbb{N}^*$ la taille de la population et $n = 0, 1, 2, \dots$ les générations successives. Sur le locus étudié, on peut trouver deux allèles différents notés $A$ et $B$. La variable aléatoire $X_n$ compte le nombre d'allèles $A$ à la génération $n$. La population à la génération $n+1$ est déduite de celle de la génération $n$ par un tirage avec remise de $N$ individus de la génération $n$ où la probabilité d'obtenir $A$ est $X_n/N$. On peut réaliser $(X_n)_{n \ge 0}$ comme la suite récurrente aléatoire définie par la relation récurrence

$$ X_{n+1} = \sum_{k=1}^{N} 1_{\{U_{n+1,k} \le \psi_{X_n}\}} \quad \text{où} \quad \psi_x := \frac{x}{N} $$

et où $(U_{n,k})_{n \ge 1, 1 \le k \le N}$ sont des variables aléatoires i.i.d. de loi uniforme sur $[0, 1]$, indépendantes de $X_0$. Pour tous $x_0, \dots, x_n$, on a

$$ \text{Loi}(X_{n+1} | X_0 = x_0, \dots, X_n = x_n) = \text{Loi}(X_{n+1} | X_n = x_n) = \text{Bin}(N, \psi_{x_n}). $$
**Remarque 12.4** (Modèle à plus de deux allèles). Le modèle de Wright-Fisher à $\ell \ge 2$ allèles $A_1, \dots, A_\ell$ fait intervenir la loi multinomiale. Plus précisément, $X_n$ prend ses valeurs dans $\{x \in \mathbb{N}^\ell : x_1 + \dots + x_\ell = N\}$ et

$$ \text{Loi}(X_{n+1} | X_0, \dots, X_n) = \text{Loi}(X_{n+1} | X_n) = \text{Mul}(N, \frac{X_n}{N}). $$

La suite $(X_n)_{n \ge 0}$ est une chaîne de Markov d'espace d'états $\{0, 1, \dots, N\}$ et de matrice de transition $\mathbf{P}$ donnée pour tous $x, y \in \{0, 1, \dots, N\}$ par

$$ \mathbf{P}(x, y) := \mathbb{P}(X_{n+1} = y | X_n = x) = \binom{N}{y} \psi_x^y (1 - \psi_x)^{N-y}. $$

Les états 0 et $N$ sont absorbants (deux classes de récurrence singleton) tandis que tous les autres états mènent à $\{0, N\}$ et sont donc transitoires. Tout comme pour le modèle de Moran, lorsque la population ne contient plus qu'un allèle, on dit qu'il est fixé dans la population. Le temps de fixation est

$$ T := \inf\{n \ge 0 : X_n \in \{0, N\}\}. $$

L'espace d'états étant fini, le temps d'atteinte $T$ de $\{0, N\}$ vérifie

$$ \mathbb{P}(T < \infty) = 1 \quad \text{et} \quad \mathbb{E}(T) < \infty. $$

De plus $(X_n)_{n \ge 0}$ converge p.s. vers la variable aléatoire $X_T$ qui suit une loi de Bernoulli portée par 0 et $N$. On a $\mathbb{P}(X_T = 0) + \mathbb{P}(X_T = N) = 1$. L'événement $\{X_T = N\}$ signifie que l'allèle $A$ est fixé tandis que l'événement $\{X_T = 0\}$ signifie que l'allèle $B$ est fixé.

$$ \text{On adopte les notations } \mathbb{P}_x = \mathbb{P}(\cdot | X_0 = x) \text{ et } \mathbb{E}_x = \mathbb{E}(\cdot | X_0 = x). $$

**Théorème 12.5** (Probabilité de fixation). Pour tout $x \in \{0, 1, \dots, N\}$,

$$ \mathbb{P}_x(X_T = N) = 1 - \mathbb{P}_x(X_T = 0) = \frac{x}{N}. $$

*Démonstration.* La preuve est identique à celle faite pour le modèle de Moran (théorème 12.2) car $\text{Loi}(X_n | X_{n-1}) = \text{Bin}(N, \frac{X_{n-1}}{N})$ d'où

$$ \mathbb{E}(X_n) = \mathbb{E}(\mathbb{E}(X_n | X_{n-1})) = \mathbb{E}(N \frac{X_{n-1}}{N}) = \mathbb{E}(X_{n-1}) = \dots = \mathbb{E}(X_0) = x, $$

□

**Remarque 12.6** (Martingale). L'intégralité de la remarque 12.3 (modèle de Moran) reste valable pour le modèle de Wright-Fisher, car pour tout $n \ge 1$,

$$ \mathbb{E}(X_n | \mathcal{F}_{n-1}) = \mathbb{E}(X_n | X_{n-1}) = N \frac{X_{n-1}}{N} = X_{n-1}, $$

où $\mathcal{F}_n := \sigma(X_0, \dots, X_n)$ pour tout $n \ge 0$.
<figure><img src="image_37.png" /><figcaption>Fig. 12.1. Quelques trajectoires, jusqu'au temps $t = 200$, de la chaîne de Wright-Fisher pour une population de taille $N = 100$, toutes issues de $X_0 = N/4 = 25$.</figcaption></figure>

**Remarque 12.7** (Couplage avec un jeu de pile ou face). Il est possible d'obtenir une minoration de $T$ par une variable aléatoire géométrique par couplage. On commence par observer que pour tout $x$ et tout $n$ on a

$$P(X_{n+1} \in \{0, N\} | X_n = x) = P(x, \{0, N\}) \ge p_*$$

où

$$p_* = \min_{0 \le x \le N} P(x, \{0, N\}) = \min_{0 \le x \le N} (1 - \psi_x)^N + \psi_x^N \ge 2^{-N+1}$$

où le minimum est atteint en $\psi_x = 1/2$. Ainsi, la chaîne atteint un état absorbant 0 ou $N$ avant qu'un lanceur de pièce n'obtienne son premier pile avec une pièce qui fait pile avec probabilité $2^{-N+1}$. Pour rendre rigoureux ce raisonnement, on construit le couple de processus (jeu de pile ou face et chaîne de Wright-Fisher) sur le même espace de probabilité. Représentons la chaîne de Markov comme une suite récurrente aléatoire : pour tout $x \in \{0, \dots, N\}$, soit $(I_{x,y})_{0 \le y \le N}$ une partition de l'intervalle $[0, 1]$, avec $I_{x,y}$ de longueur

$$|I_{x,y}| = P(x, y),$$

et soit $(U_n)_{n \ge 1}$ une suite de variables aléatoires i.i.d. de loi uniforme sur $[0, 1]$, indépendante de $X_0$. On réalise $(X_n)_{n \ge 0}$ en utilisant la récurrence

$$X_{n+1} := y \quad \text{où } y = y(X_n, U_{n+1}) \text{ est tel que } U_{n+1} \in I_{X_n, y}.$$
Il est possible de faire un dessin figurant une tour de $N+1$ copies de l'intervalle $[0,1]$ avec ses $N+1$ partitions. Dans cette construction, on peut toujours ordonner les intervalles, et supposer que $I_{x,0}$ et $I_{x,N}$ sont les deux premiers en partant de la gauche. Par construction $p_* = \min_{0 \le x \le N} |I_{x,0} \cup I_{x,N}|$ et

$$T := \inf\{n \ge 0 : X_n \in \{0, N\}\} \le T_* := \inf\{n \ge 1 : 1_{\{U_n \le p_*\}} = 1\}$$

et $T_*$ suit la loi géométrique $\text{Geo}(p_*)$. En particulier, $T$ est intégrable et donc fini p.s. Cette méthode de couplage fournit plus généralement, pour des chaînes de Markov, un minorant géométrique des temps d'atteinte d'ensembles.

La vitesse de fixation peut être quantifiée par l'hétérozygotie, c'est-à-dire la probabilité que deux gènes choisis aléatoirement et sans remise dans la population totale à la génération $n$ soient représentés par des allèles différents. Elle est donnée par

$$H_n = \frac{2X_n(N - X_n)}{N(N - 1)} = \frac{2}{N - 1} \text{Var}(X_{n+1} | X_n).$$

On a $H_n = 0$ si $n \ge T$ et donc $H_n \xrightarrow{n \to \infty} 0$. Le théorème suivant montre que l'hétérozygotie moyenne $\mathbb{E}(H_n)$ décroît exponentiellement au cours du temps.

**Théorème 12.8 (Hétérozygotie moyenne).** Pour tout $n \ge 1$,

$$\mathbb{E}(H_n) = h_0\lambda^n \quad \text{où } h_0 := \mathbb{E}(H_0), \quad \lambda := 1 - \frac{1}{N} \in ]0, 1[,$$

et

$$\text{Var}(X_n) = \mathbb{E}(X_0)(N - \mathbb{E}(X_0))(1 - \lambda^n) + \lambda^n \text{Var}(X_0).$$

**Démonstration.** Il suffit d'établir que pour tout $n \ge 1$,

$$\mathbb{E}(X_n(N - X_n)) = \left(1 - \frac{1}{N}\right)^n \mathbb{E}(X_0(N - X_0)).$$

On a

$$\mathbb{E}(X_n(N - X_n)) = N\mathbb{E}(X_n) - \mathbb{E}(X_n^2) = N\mathbb{E}(X_{n-1}) - \mathbb{E}(\mathbb{E}(X_n^2 | X_{n-1})).$$

On écrit alors

$$\begin{aligned} \mathbb{E}(X_n^2 | X_{n-1}) &= \text{Var}(X_n | X_{n-1}) + (\mathbb{E}(X_n | X_{n-1}))^2 \\ &= X_{n-1} \left(1 - \frac{X_{n-1}}{N}\right) + X_{n-1}^2. \end{aligned}$$

En regroupant les termes, on obtient bien

$$\mathbb{E}(X_n(N - X_n)) = \left(1 - \frac{1}{N}\right) \mathbb{E}(X_{n-1}(N - X_{n-1})).$$
ce qui fournit la formule pour $E(H_n)$. Pour la variance, on utilise la formule<fnref n="4" />

$$ \mathrm{Var}(X_n) = \mathrm{E}(\mathrm{Var}(X_n | X_{n-1})) + \mathrm{Var}(\mathrm{E}(X_n | X_{n-1})). $$

□

Comme dans le problème de la ruine du joueur (théorème 2.2), et comme pour le modèle de Moran (théorème 12.2), on souhaite déterminer le temps moyen de sortie en fonction du point de départ. Pour $t = 0 \le x \le N$, notons $m(x) = E_x(T)$. Alors $m(0) = m(N) = 0$ et, grâce à la propriété de Markov (faible), on obtient un système d'équations linéaires : pour tout $0 < x < N$,

$$ m(x) = \sum_{y=0}^{N} E_x(T1_{\{X_1=y\}}) = \sum_{y=0}^{N} (1 + E_y(T))\mathbf{P}(x, y) = 1 + \sum_{y=0}^{N} m(y)\mathbf{P}(x, y), $$

dont $m$ est l'unique solution positive minimale. Bien qu'on puisse calculer $m$ numériquement, il n'existe pas d'expression simple et explicite de $m$ comme pour le temps de sortie d'un intervalle pour la marche aléatoire simple (théorème 2.2) ou pour le modèle de Moran (théorème 12.2). Il est toutefois possible de trouver l'équivalent de $m$ lorsque la taille de la population $N$ tend vers l'infini : pour tout $0 < p < 1$,

$$ m([Np]) \sim_{N \to \infty} -2N(p \log(p) + (1-p) \log(1-p)), $$

ce qui rappelle le modèle de Moran (théorème 12.2). La clé de la démonstration de ce résultat consiste à établir que la suite de processus $Y^{(N)}$, définie par

$$ Y^{(N)} = \left( \frac{1}{N} X_{\lfloor Nt \rfloor} \right)_{t \in \mathbb{R}_+} $$

à valeurs dans $[0, 1]$, converge en loi quand $N \to \infty$ vers un processus de diffusion $Y$ sur $[0, 1]$, solution de l'équation différentielle stochastique

$$ dY_t = \sqrt{Y_t(1 - Y_t)} dB_t. $$

On trouvera les détails de ce résultat dans le chapitre 27.

**Remarque 12.9** (Lien avec le modèle de Moran). Pour le modèle de Wright-Fisher, l'espérance du temps d'absorption est de l'ordre de $N$ tandis qu'il est de l'ordre de $N^2$ pour le modèle de Moran (voir théorème 12.2). Cela vient du fait que lors d'une transition du modèle de Moran, un seul allèle est modifié tandis que tous sont concernés à chaque transition du modèle de Wright-Fisher.

<footnote n="4">4. Plus généralement, si $\Phi$ est convexe alors $\mathrm{Var}^{\Phi}(X) := \mathrm{E}(\Phi(X)) - \Phi(\mathrm{E}(X)) \ge 0$ par l'inégalité de Jensen, et $\mathrm{Var}^{\Phi}(X) = \mathrm{E}(\mathrm{Var}^{\Phi}(X|Y)) + \mathrm{Var}^{\Phi}(\mathrm{E}(X|Y))$ où $\mathrm{Var}^{\Phi}(X|Y) := \mathrm{E}(\Phi(X)|Y) - \Phi(\mathrm{E}(X|Y))$. Si $\Phi(t) = t^2$ on retrouve la variance.</footnote>
## 12.4 Modèle de Wright-Fisher avec mutations

Supposons à présent qu'avant le tirage avec remise qui permet de passer d'une génération à l'autre, chaque individu puisse muter indépendamment des autres : l'allèle A mute en allèle B avec probabilité $u \in [0, 1]$ et B mute en A avec probabilité $v \in [0, 1]$. Le passage d'une génération à la suivante se fait donc en quelque sorte en appliquant d'abord un opérateur de mutation, puis un opérateur d'héritage (procéder dans l'ordre inverse conduirait à un processus différent). Ce mécanisme définit un nouveau processus $(X_n)_{n \ge 0}$ à valeurs dans $\{0, 1, \dots, N\}$ qui vérifie

$$ \text{Loi}(X_{n+1} | X_0, \dots, X_n) = \text{Bin}(N, \psi_{X_n}) $$

où à présent

$$ \psi_x = \frac{x(1-u) + (N-x)v}{N}. $$

Le processus X est une chaîne de Markov d'espace d'états $\{0, 1, \dots, N\}$.

- si $u = v = 0$, alors on retrouve le modèle sans mutation étudié précédemment. Dans ce cas, les états 0 et $N$ sont absorbants tandis que $\{1, \dots, N-1\}$ est une classe transitoire, l'ensemble des lois invariantes est $\{(1-p)\delta_0 + p\delta_N : p \in [0, 1]\}$, et partant de $x$, la chaîne converge en loi vers $(1-p_x)\delta_0 + p_x\delta_N$ où $p_x = x/N$;

- si $u = 0$ et $0 < v \le 1$ alors l'état $N$ est absorbant tandis que l'ensemble d'états $\{0, \dots, N-1\}$ est une classe transitoire, et la chaîne converge en loi vers l'unique loi invariante $\delta_N$;

- si $0 < u \le 1$ et $v = 0$ alors l'état 0 est absorbant tandis que l'ensemble d'états $\{1, \dots, N\}$ est une classe transitoire, et la chaîne converge en loi vers l'unique loi invariante $\delta_0$;

- si $u = v = 1$, alors $\{0, N\}$ est une classe de récurrence de période 2 tandis que l'ensemble d'états $\{1, \dots, N-1\}$ est une classe transitoire. La chaîne est presque sûrement absorbée par la classe $\{0, N\}$ et oscille ensuite périodiquement entre les états 0 et $N$;

- si $0 < u \le 1$ et $0 < v < 1$ ou si $0 < u < 1$ et $0 < v \le 1$ alors la chaîne est irréductible. Comme l'espace d'états est fini, elle est récurrente positive et possède une unique loi invariante $\mu$, et cette loi charge tous les états. La loi des grands nombres pour les chaînes de Markov indique que quelle que soit la loi initiale, p.s. pour tout $x \in \{0, 1, \dots, N\}$, on a

$$ \lim_{n \to \infty} \frac{\text{card}\{0 \le k \le n : X_k = x\}}{n+1} = \mu(x). $$

Comme de plus la matrice de transition possède un coefficient diagonal non nul, la chaîne est apériodique, et converge donc en loi vers $\mu$, quelle que soit la loi initiale, c'est-à-dire que pour tout $x \in \{0, 1, \dots, N\}$,

$$ \lim_{n \to \infty} \mathbb{P}(X_n = x) = \mu(x). $$
La loi $\mu$ est calculable numériquement, mais ne possède pas de formulation simple. Cependant, il est possible de déterminer explicitement sa moyenne $m$ par exemple. En effet, comme $\mu$ est invariante,

$$m = \sum_{x=0}^{N} y\mu(y) = \sum_{y=0}^{N} y \sum_{x=0}^{N} \mu(x) \mathbf{P}(x, y) = \sum_{x=0}^{N} \mu(x) \sum_{y=0}^{N} y \mathbf{P}(x, y).$$

Or comme $\mathbf{P}(x, \cdot) = \text{Bin}(N, \frac{(1-u)x+v(N-x)}{N})$, on a

$$\sum_{y=0}^{N} y \mathbf{P}(x, y) = (1-u)x + v(N-x).$$

Cette formule, affine en $x$, donne une équation fermée pour $m$ :

$$m = \sum_{x=0}^{N} ((1-u)x + v(N-x))\mu(x) = (1-u)m + v(N-m),$$

d'où enfin

$$m = \frac{Nv}{u+v}.$$

Pour la variance $\sigma^2$ de $\mu$, un calcul permet d'établir que

$$\sigma^2 = \frac{N^2 uv}{(u+v)^2(2N(u+v)+1)} + o_{N \to \infty}(N).$$

Notons que si $u = v$ alors $m = N/2$ ne dépend plus de $u$ et $v$. Au-delà des deux premiers moments, on peut établir, comme présenté dans le chapitre 27, que le processus

$$Y^{(N)} = \left( \frac{1}{N} X_{\lfloor Nt \rfloor} \right)_{t \in \mathbb{R}_+}$$

converge en loi quand $(N, Nu, Nv) \to (\infty, \alpha, \beta)$ vers un processus de diffusion $Y$ sur $[0, 1]$, solution de l'équation différentielle stochastique

$$dY_t = \sqrt{Y_t(1-Y_t)}dB_t - \alpha Y_t dt + \beta (1-Y_t)dt,$$

dont la loi invariante est une loi Beta$(2\beta, 2\alpha)$, dont la moyenne et la variance sont données par

$$\frac{\beta}{\alpha + \beta} \quad \text{et} \quad \frac{\alpha\beta}{(\alpha + \beta)^2(2(\alpha + \beta) + 1)}.$$
<figure><img src="image_38.png" /><figcaption>**Fig. 12.2.** Loi invariante et histogramme d'un échantillon de taille 10000 de $X_{100}$ avec $X_0 = 50$, pour une population de taille $N = 100$, avec taux de mutation $u = v = 0.05$.</figcaption></figure>

## 12.5 Modèle de Wright-Fisher avec sélection

Sur le plan biologique, il est possible que les allèles A et B ne procurent pas à l'individu qui en est doté des capacités identiques : viabilité, potentiel attractif, etc. Cette asymétrie est responsable de la sélection naturelle. Considerons pour simplifier le modèle sans mutation. Il est possible de prendre en compte un avantage sélectif de l'allèle A en effectuant la transformation suivante de la génération n avant la fabrication de la génération n + 1 : chaque individu de type A est remplacé par (1 + s) individus de type A, où s > 0 est un paramètre fixé (penser à une urne). Cela revient à prendre

$$ \text{Loi}(X_{n+1} | X_0, \dots, X_n) = \text{Loi}(X_{n+1} | X_n) = \text{Bin}(N, \psi_{X_n}) $$

avec cette fois-ci

$$ \psi_x = \frac{x(1 + s)}{x(1 + s) + N - x}. $$

Le processus $(X_n)_{n \ge 0}$ est une chaîne de Markov. Son espérance ne vérifie pas une équation de récurrence linéaire car $\psi_x$ n'est pas linéaire en $x$ (contrairement au cas sans sélection $s = 0$). Les états 0 et $N$ sont absorbants, tandis que $\{1, \dots, N-1\}$ est une classe transitoire. Comme $\psi_x$ n'est plus linéaire en $x$, on ne peut plus calculer la probabilité de fixation avec la méthode utilisée dans le cas sans sélection ($s = 0$). Soit $\mathbf{P}$ la matrice de transition de $X$, et

$$ \pi(x) := \mathbf{P}_x(X_T = N) $$
la probabilité de fixation de l'allèle A. On a

$$ \pi(x) = \sum_{y=0}^{N} \mathbf{P}(x, y)\pi(y). $$

Il est possible de trouver un équivalent de $\pi$ lorsque la taille de la population $N$ tend vers l'infini et $s$ est proche de 0. Plus précisément, on montre dans le chapitre 27 que si $Ns \to \alpha > 0$ alors, pour tout $0 < p < 1$,

$$ \pi(\lfloor pN \rfloor) \underset{N \to \infty}{\longrightarrow} \frac{1 - e^{-2\alpha p}}{1 - e^{-2\alpha}}. $$

Ainsi, le passage de $s = 0$ à $s > 0$ modifie considérablement les probabilités de fixation.

## 12.6 Modèle de Cannings

Dans le modèle de Wright-Fisher, on passe de la génération $n$ à la génération $n+1$ en effectuant un tirage avec remise de $N$ individus dans la génération $n$. Cela revient à attribuer un certain nombre d'enfants à chaque individu de la génération $n$, les enfants héritant de l'allèle du père. Si $y_1, \dots, y_N$ désignent le nombre d'enfants de chaque individu de la génération $n$ (numérotés de 1 à $N$) alors on a $y_1 + \dots + y_N = N$ et $(y_1, \dots, y_N)$ suit la loi multinomiale symétrique de taille $N$ et de paramètre $(1/N, \dots, 1/N)$. Cela correspond à lancer $N$ fois un dé à $N$ faces équilibrées. Ce mécanisme qui associe chaque individu de la génération $n+1$ à un individu de la génération $n$ peut être utilisé pour étudier la transmission de gènes quel que soit le nombre d'allèles considéré : les enfants héritent de l'allèle du père. Le modèle de Cannings va au-delà de cette structure multinomiale symétrique en permettant une dépendance entre les individus d'une même génération : on suppose seulement que le vecteur $(y_1, \dots, y_N)$ vérifie $y_1 + \dots + y_N = N$ et possède une loi échangeable, c'est-à-dire que pour toute permutation $\tau \in S_N$,

$$ (y_1, \dots, y_N) \stackrel{\text{loi}}{=} (y_{\tau(1)}, \dots, y_{\tau(N)}). $$

En particulier, les $y_1, \dots, y_N$ ont même loi. Comme $y_1 + \dots + y_N = N$ on en déduit que $E(y_1) = 1$ et on note $\sigma^2 = \text{Var}(y_1)$. On en déduit également $\text{Cov}(y_1, y_2) = -\sigma^2/(N-1)$ car

$$ 0 = \text{Var}(y_1 + \dots + y_N) = N\text{Var}(y_1) + N(N-1)\text{Cov}(y_1, y_2). $$

**Remarque 12.10** (Modèle de Moran). Tout comme le modèle de Wright-Fisher, le modèle de Moran est un cas particulier du modèle de Cannings. Il correspond à la loi de $(y_1, \dots, y_N)$ suivante : la transition de la génération $n$ à la génération $n+1$ se fait en tirant avec remise deux individus dans la génération $n$, le premier a deux fils, le second aucun, et tous les autres individus ont exactement un seul descendant.
Dans le cas à deux allèles, ce mécanisme de transition donne lieu, tout comme pour le modèle de Wright-Fisher, à une chaîne de Markov $(X_n)_{n\ge0}$ d'espace d'états $\{0,1,\dots,N\}$ donnée pour tout $x \in \{0,1,\dots,N\}$ par

$$ \mathrm{Loi}(X_{n+1} | X_n = x) = \mathrm{Loi}(y_1 + \dots + y_x). $$

Les états 0 et $N$ sont absorbants tandis que tous les autres mènent à $\{0, N\}$ et sont donc transitoires. Notons que $\mathrm{E}(X_{n+1} | X_n = x) = x : (X_n)_{n\ge0}$ est donc une martingale. De plus,

$$ \begin{aligned} \mathrm{Var}(X_{n+1} | X_n = x) &= \mathrm{Var}(y_1 + \cdots + y_x) \\ &= x\sigma^2 + x(x-1)\mathrm{Cov}(y_1, y_2) \\ &= x\sigma^2 - \frac{x(x-1)\sigma^2}{N-1} \\ &= \sigma^2 \frac{x(N-x)}{N-1}. \end{aligned} $$

**Remarque 12.11** (Loi multinomiale). Soit $n, d \in \mathbb{N}^*$ et $(p_1, \dots, p_d) \in [0, 1]^d$ avec $p_1 + \cdots + p_d = 1$. La loi multinomiale de taille $n$ et de paramètre $(p_1, \dots, p_d)$ est la loi sur l'ensemble fini

$$ \{(n_1, \dots, n_d) \in \mathbb{N}^d : n_1 + \cdots + n_d = n\} $$

donnée par $(e_1, \dots, e_d)$ est la base canonique de $\mathbb{R}^d$)

$$ (p_1\delta_{e_1} + \cdots + p_d\delta_{e_d})^{*n} = \sum_{\substack{(n_1, \dots, n_d) \in \mathbb{N}^d \\ n_1 + \cdots + n_d = n}} \frac{n!}{n_1! \cdots n_d!} p_1^{n_1} \cdots p_d^{n_d} \delta_{(n_1, \dots, n_d)}. $$

Elle code n lancers d'un dé à d faces avec probabilités d'apparition $p_1, \dots, p_d$. Elle est «symétrique» lorsque $p_1 = \cdots = p_d = 1/d$. Si $(X_1, \dots, X_d)$ suit la loi $\mathrm{Mul}(n, (p_1, \dots, p_d))$ et si $I_1, \dots, I_r$ est une partition de $\{1, \dots, n\}$ alors

$$ (X_{I_1}, \dots, X_{I_r}) \sim \mathrm{Mul}(n, (p_{I_1}, \dots, p_{I_r})) $$

où

$$ X_I := \sum_{i \in I} X_i \quad \text{et} \quad p_I := \sum_{i \in I} p_i. $$

Si $I \subset \{1, \dots, n\}$ alors $X_I \sim \mathrm{Bin}(n, p_I)$. Ainsi les lois marginales d'une loi multinomiale sont binomiales. On identifie la loi $\mathrm{Bin}(n, p)$ à $\mathrm{Mul}(n, (p, 1-p))$. Si $(X_1, \dots, X_d) \sim \mathrm{Mul}(n, (p_1, \dots, p_d))$ alors $X_d = n - X_{\{1, \dots, d-1\}}$.

**Théorème 12.12** (Spectre de la matrice de transition). Les valeurs propres de la matrice de transition du modèle de Cannings sont données par

$$ \lambda_0 = 1 \quad \text{et} \quad \lambda_j = \mathrm{E}(y_1 y_2 \cdots y_j) \quad \text{pour} \quad 1 \le j \le N. $$

En particulier, pour le modèle de Wright-Fisher, $\lambda_0 = 1$ et, pour $1 \le j \le N$,

$$ \lambda_j = \frac{N(N-1)\cdots(N-j+1)}{N^j} = \left(1 - \frac{1}{N}\right) \cdots \left(1 - \frac{j-1}{N}\right). $$
Démonstration. Considérons la matrice de Vandermonde

$$
V = \begin{pmatrix}
1 & 0 & 0 & \cdots & 0 \\
1 & 1 & 1^2 & 1^3 & \cdots & 1^N \\
1 & 2 & 2^2 & 2^3 & \cdots & 2^N \\
1 & 3 & 3^2 & 3^3 & \cdots & 3^N \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
1 & N & N^2 & N^3 & \cdots & N^N
\end{pmatrix} \in \mathcal{M}_{N+1,N+1}(\mathbb{R}).
$$

Soit $\mathbf{P} \in \mathcal{M}_{N+1,N+1}(\mathbb{R})$ la matrice de transition du modèle de Cannings.
Pour tous $i,j \in \{0,1,\dots,N\}$, on a

$$
(\mathbf{P}V)_{i,j} = \sum_{k=0}^{N} \mathbf{P}_{i,k} V_{k,j} = \sum_{k=0}^{N} \mathbf{P}_{i,k} k^j = \mathbb{E}(X_1^j | X_0 = i).
$$

Par ailleurs, comme

$$
\mathrm{Loi}(X_1 | X_0 = i) = \mathrm{Loi}(y_1 + \cdots + y_i),
$$

il existe des réels $b_{i,k}$ tels que

$$
\mathrm{E}(X_1^j | X_0 = i) = \mathrm{E}((y_1 + \dots + y_i)^j) = i^{[j]} \mathrm{E}(y_1 y_2 \dots y_j) + \sum_{k=0}^{j-1} b_{i,k} i^k,
$$

où

$$
i^{[j]} := i(i-1)(i-2)\cdots(i-j+1).
$$

Pour le voir, on peut procéder par récurrence sur $j$ en observant que la loi $\text{Loi}(y_1, \dots, y_i | y_k)$ est échangeable. Il existe donc une matrice triangulaire supérieure $T \in \mathcal{M}_{N+1,N+1}(\mathbb{R})$ telle que $T_{0,0} = 1$, $T_{j,j} = \text{E}(y_1 \cdots y_j)$ pour tout $j \in \{1, \dots, N\}$, et pour tous $i, j \in \{0, 1, \dots, N\}$,

$$
\mathrm{E}(X_1^j | X_0 = i) = \sum_{k=0}^{j} T_{k,j} i^k = (VT)_{i,j}.
$$

On a donc

$$
\mathbf{P}V = VT,
$$

et comme $V$ est inversible, les matrices $\mathbf{P}$ et $T$ ont même spectre.

<figure><img src="image_39.png" /></figure>

On a $\mathrm{E}(y_1) = 1$ et donc $1$ est valeur propre de multiplicité $2$, ce qui correspond au fait que la chaîne possède deux classes de récurrence $\{0\}$ et $\{N\}$. La deuxième plus grande valeur propre est

$$
\lambda_2 = \mathrm{E}(y_1 y_2) = \mathrm{Cov}(y_1, y_2) + \mathrm{E}(y_1)^2 = 1 - \frac{\sigma^2}{N-1}.
$$
12.7 Pour aller plus loin

Les mécanismes de la génétique ont été considérablement discutés depuis leur découverte par Gregor Mendel au dix-neuvième siècle. La loi de Godfrey Harold Hardy et de Wilhelm Weinberg a été découverte de manière indépendante par ces deux scientifiques au tout début du vingtième siècle. Le modèle de Ronald Aylmer Fisher et de Sewall Wright date des années 1930. Le modèle de John Moran date des années 1950 tandis que le modèle de Chris Cannings date des années 1970. Ces modèles et leurs extensions sont étudiés par exemple dans les livres de Jean-François Delmas et Benjamin Jourdain [DJ06], de Rick Durrett [Dur08], de Warren Ewens [Ewe04], de Simon Tavaré [Tav04], ainsi que dans le cours de Sylvie Méléard [Mé13]. Ils inspirent encore aujourd'hui des recherches mêlant biologie et mathématiques. L'approximation des chaînes de Markov par des diffusions fait l'objet du chapitre 27. La généalogie du processus de Wright-Fisher est étudiée dans le chapitre 13, qui mène au processus des restaurants chinois du chapitre 14.
# Généalogies et coalescence

**Mots-clés.** Arbre; branchement; généalogie; coalescent de Kingman.

**Outils.** Loi géométrique; loi exponentielle; loi multinomiale; loi de Poisson; loi d'Ewens; processus de Galton-Watson; modèle de Wright-Fisher.

**Difficulté.** **

Dans ce chapitre, on étudie la généalogie du processus de Wright-Fisher du chapitre 12. On la relie notamment au processus des restaurants chinois du chapitre 14. Cette étude conduit au processus de coalescence de Kingman.

On s'intéresse à l'évolution d'une population de taille finie $N$ constante au fil des générations. On suppose que les générations ne se chevauchent pas, et que pour tout $n \ge 0$, la génération $n$ meurt en donnant naissance à la génération $n+1$. Le mécanisme de transition est markovien : la génération $n+1$ dépend de la génération $n$ et d'une source d'aléa indépendante. À chaque génération, on numérote les individus de 1 à $N$. Pour tout $n \ge 0$, soit $a_i^{n+1}$ le numéro du parent de l'individu $i$ de la génération $n+1$. Ce parent appartient à la génération $n$. On suppose que les v.a. $a_1^{n+1}, \dots, a_N^{n+1}$ sont i.i.d. de loi uniforme sur $\{1, \dots, N\}$. On note $\nu_i^n$ le nombre de descendants à la génération $n+1$ de l'individu $i$ vivant à la génération $n$. Les v.a. $(\nu_i^n)_{1 \le i \le N}$ ne sont pas indépendantes puisqu'elles vérifient la relation $\nu_1^n + \dots + \nu_N^n = N$. On a

$$ \nu_i^n = \sum_{j=1}^{N} 1_{\{a_j^{n+1}=i\}}. $$

Par conséquent, le vecteur aléatoire $\nu^n = (\nu_1^n, \dots, \nu_N^n)$ suit la loi multinomiale de taille $N$ et de paramètre $(1/N, \dots, 1/N)$, c'est-à-dire

$$ \mathbb{P}(\nu_1^n = n_1, \dots, \nu_N^n = n_N) = \frac{N!}{n_1! \dots n_N!} \left(\frac{1}{N}\right)^N 1_{\{n_1+\dots+n_N=N\}}. $$
Si chaque individu est réduit à un allèle A ou B d'un gène à deux allèles et si le nombre d'allèles A à la génération $n$ est égal à $x$, alors la loi du nombre d'allèles A à la génération $n+1$ suit une loi binomiale $Bin(N, x/N)$. On retrouve donc le modèle de Wright-Fisher du chapitre 12. Le lien avec la loi multinomiale a déjà été observé lors de l'étude du modèle de Cannings dans le chapitre 12. Nous avons construit ici le modèle de Wright-Fisher à rebours en remontant les générations, comme illustré dans la figure 13.1. Chaque individu de la génération $n+1$ possède un ancêtre à la génération $n$ choisi indépendamment des autres et selon la loi uniforme sur $\{1, ..., N\}$. Certains membres de la génération $n$ n'ont donc aucun descendant à la génération $n+1$ et c'est ce qui va entraîner le phénomène dit de dérive génétique (fixation d'un allèle).

**Remarque 13.1** (Lois de Poisson et multinomiale). Si $Y_1, \dots, Y_N$ sont des v.a. indépendantes de loi de Poisson de paramètres respectifs $\theta_1, \dots, \theta_N$ alors

$$ \text{Loi}((Y_1, \dots, Y_N) | Y_1 + \dots + Y_N = N) = \text{Mul}(N, (p_1, \dots, p_N)) $$

où $p_k = \theta_k/(\theta_1 + \dots + \theta_N)$ pour tout $1 \le k \le N$. En particulier, dans le cas où $\theta_1 = \dots = \theta_N$ alors $p_1 = \dots = p_N = 1/N$. Ainsi le modèle de Wright-Fisher revient à un processus de Galton-Watson de loi de reproduction Poisson, conditionné à avoir une taille de population fixée, voir chapitre 3.

Le modèle de Wright-Fisher sans mutation est connu pour faire apparaître un phénomène dérive génétique, c'est-à-dire qu'un allèle l'emporte sur l'autre et que la population devient homozygote après un nombre aléatoire mais fini presque sûrement de générations. De même, après un certain nombre de générations, tous les individus proviendront donc d'un même ancêtre. C'est ce phénomène qu'il s'agit à présent de décrire et quantifier.

Comme les individus donnés choisissent leur parent de manière i.i.d. et uniforme, la probabilité que deux individus fixés aient deux parents distincts (à la génération précédente) est égale à $N(N-1)/N^2 = 1-1/N$. Ces parents choisissant eux-mêmes leurs parents indépendamment, la probabilité pour que le premier ancêtre commun à deux individus donnés remonte à plus de $r$ générations vaut $(1-1/N)^r$. Le temps $T_2$ d'apparition de l'ancêtre commun le plus récent (ACPR) suit donc une loi géométrique sur $\mathbb{N}^*$ de paramètre $p_2 = 1/N$. Si l'on considère à présent trois individus, notons $T'_3$ l'instant d'apparition de l'ACPR pour deux individus au moins (premier temps de coalescence parmi trois individus). On a

$$ P(T'_3 > 1) = \frac{N(N-1)(N-2)}{N^3} = 1 - \frac{3}{N} + \frac{2}{N^2}. $$

Les choix des parents sont indépendants à chaque génération, d'où, pour $k \ge 1$,

$$ P(T'_3 > k) = \left(1 - \frac{3}{N} + \frac{2}{N^2}\right)^k. $$

La loi de $T'_3$ est donc la loi géométrique de paramètre $p_3 = 3/N - 2/N^2$. On note $T_3 = T'_3 + T'_2$ le temps d'apparition d'un ACPR pour les trois individus
<figure><img src="image_40.png" /><figcaption>**Fig. 13.1.** Dans le modèle de Wright-Fisher, tout se passe comme si pour chaque génération, chaque individu, indépendemment de tous les autres, choisissait uniformément un père dans la génération précédente et en héritait. Le graphique ci-dessus illustre l'évolution de ces relations de filiation sur quelques générations pour une population de taille $N = 5$. Ce point de vue généalogique est au cœur de ce chapitre.</figcaption></figure>

(ici $T'_2 := T_3 - T'_3$). Remarquons que si l'ACPR de deux individus est l'ACPR des trois individus alors $T'_2 = 0$. Déterminons la loi de $T'_2$. La probabilité que trois individus distincts aient le même père est $1/N^2$. La probabilité que $T'_2$ soit nul sachant que $T'_3 = k$ est donc égale à

$$ \mathbb{P}(T'_{2} = 0 | T'_{3} = k) = \frac{\mathbb{P}(T'_{2} = 0, T'_{3} = k)}{\mathbb{P}(T'_{3} = k)} = \frac{(1 - p_{3})^{k-1}(1/N^{2})}{p_{3}(1 - p_{3})^{k-1}} = \frac{1}{3N - 2}. $$

Comme les choix des parents sont indépendants d'une génération à l'autre, on déduit de la première partie que la loi de $T'_2$ sachant que $T'_3 = k$ et $T'_2 > 0$ est la loi géométrique de paramètre $p_2$. Cela implique que sachant $T'_3 = k$, la variable $T'_2$ est égale en loi à $UV$ où $U$ et $V$ sont indépendantes de lois respectives $Ber(1 - 1/(3N - 2))$ et $Geo_{N^*}(p_2)$.

Plus généralement, le nombre de parents distincts d'un groupe de $k$ individus peut être vu comme le nombre d'urnes occupées après que l'on a lancé $k$ balles, indépendamment et à chaque fois uniformément, dans $N$ urnes. Pour tout $j \in \{1, \dots, k\}$, la probabilité pour que ce nombre soit $j$ est

$$ g_{k,j}^{(N)} = \mathbb{P}(k \text{ individus ont } j \text{ parents distincts}) = \frac{N(N-1)\cdots(N-j+1)S_{k,j}}{N^k}, $$
où $S_{k,j}$ est le nombre de Stirling de seconde espèce<fnref n="1" />. En effet, il y a exactement $N(N-1)\cdots(N-j+1)$ façons de choisir $j$ parents distincts parmi $N$, et $S_{k,j}$ façons d'associer à ces $j$ parents $k$ enfants, et enfin, $N^k$ est le nombre de façons d'assigner $k$ enfants à leurs parents. On définit à présent le processus ancestral $(A_n^N(r))_{r \in \mathbb{N}}$ en notant $A_n^N(r)$ le nombre d'ancêtres distincts à la génération $-r$ pour un groupe de taille $n$ au temps 0 (le temps remonte ici). La suite $(A_n^N(r))_{r \in \mathbb{N}}$ est une chaîne de Markov sur $\{1, \dots, n\}$ de matrice de transition

$$G_N = \left( g_{j,k}^{(N)} 1_{\{j \le k\}} \right)_{1 \le j,k \le n}.$$

Pour cette chaîne, l'état 1 est absorbant tandis que les états 2, ..., $n$ sont transitoires puisqu'ils mènent tous à 1. Il est difficile d'étudier les propriétés fines de cette chaîne, comme par exemple des propriétés sur le temps d'atteinte de l'état 1. Nous allons donc remplacer ce modèle à temps discret par un modèle plus simple à temps continu.

## 13.1 Modèle généalogique à temps continu

Considérons la chaîne de Markov à temps continu $(A(t))_{t \ge 0}$ sur $\mathbb{N}^*$ dont le générateur infinitésimal est donné pour tous $j, k \in \mathbb{N}^*$ par

$$L(j,k) = \begin{cases} \binom{j}{2} & \text{si } j \ge 2 \text{ et } k = j-1, \\ -\binom{j}{2} & \text{si } j \ge 2 \text{ et } k = j, \\ 0 & \text{sinon.} \end{cases}$$

C'est un processus de mort pur sur $\mathbb{N}^*$ (sauts de $n$ à $n-1$ seulement) pour lequel 1 est absorbant. Le temps de séjour en $n$ suit la loi exponentielle de paramètre $(\frac{n}{2})$. Ceci s'interprète de la manière suivante : chaque couple d'individus se cherche un père indépendamment des autres couples et y arrive au bout d'un temps exponentiel de paramètre 1. Pour $n$ individus, on a $(\frac{n}{2})$ couples distincts. Or le minimum de $(\frac{n}{2})$ variables aléatoires indépendantes de même loi exponentielle Exp(1) suit la loi Exp$(\binom{n}{2})$ (voir aussi lemme 11.2). On note $(A_n(t))_{t \ge 0}$ le processus issu de $n$, qui est à valeurs dans $\{1, \dots, n\}$.

Montrons à présent que le processus à temps discret renormalisé converge vers le processus à temps continu. On s'intéresse au cas limite où $N$ est grand et où l'unité de temps se compte en $N$ générations. La date d'apparition de l'ACPR de deux individus donnés est donc $T'_2/N$ où $T'_2$ suit la loi géométrique sur $\mathbb{N}^*$ de paramètre $p_2 = 1/N$.

**Lemme 13.2** (Loi exponentielle comme limite de lois géométriques renormalisées). Si $(V_n)_{n \ge 1}$ est une suite de variables aléatoires de loi géométrique de paramètres respectifs $(\mu_n)_{n \ge 1}$ telle que $\lim_{n \to \infty} n\mu_n = \mu > 0$ alors $(n^{-1}V_n)_{n \ge 1}$ converge en loi vers la loi exponentielle de paramètre $\mu$.

<footnote n="1">1. Nombre de façons de découper un ensemble à $k$ éléments en $j$ ensembles non vides. Apparaît aussi dans le chapitre 1 pour étudier le collectionneur de coupons!</footnote>
Démonstration. On utilise les fonctions caractéristiques. Si $V \sim \text{Exp}(\mu)$, alors

$$ \varphi_{V_n/n}(t) = \frac{\mu_n e^{it/n}}{1 - (1 - \mu_n)e^{it/n}} = \frac{n\mu_n}{n\mu_n - n(1 - e^{it/n})} \underset{n \to \infty}{\longrightarrow} \frac{\mu}{\mu - it} = \varphi_V(t), $$

Alternativement, pour tout $x \ge 0$,

$$ \mathbb{P}(n^{-1}V_n \ge x) = \mathbb{P}(V_n \ge nx) = e^{-\mu_n \lfloor nx \rfloor} \to e^{-\mu x}. $$

□

Pour une population comportant $N$ individus et dans une échelle de temps où l'unité de temps est égale à $N$ générations, le temps $T_2$ d'apparition de l'ACPR pour deux individus donnés est approximativement de loi exponentielle de paramètre 1. En particulier $\mathbb{E}(T_2) = 1$. On a donc montré que

$$ (A_2^N(\lfloor Nt \rfloor))_{t \ge 0} \xrightarrow[N \to \infty]{\text{loi}} (A_2(t))_{t \ge 0}. $$

Considérons l'arbre généalogique de trois individus. Quand $N$ tend vers l'infini, $T_2'/N$ converge en loi vers $\text{Exp}(1)$. De plus, on montre de même que $(T_2'/N, T_3'/N)$ converge vers $(T_2, T_3)$ où $T_2$ et $T_3$ sont indépendantes de lois exponentielles respectives $\text{Exp}(1)$ et $\text{Exp}(3)$ :

$$ \begin{aligned} \varphi_{(T'_2/N, T'_3/N)}(s, t) &= \mathbb{E}(\exp(isT'_2/N + itT'_3/N)) \\ &= \mathbb{E}(\mathbb{E}(\exp(isT'_2/N)|T'_3) \exp(itT'_3/N)) \\ &= \frac{1}{3N-2} \mathbb{E}(e^{itT'_3/N}) + \frac{3N-3}{3N-2} \mathbb{E}(e^{itV/N}) \mathbb{E}(e^{itT'_3/N}) \end{aligned} $$

et le lemme 13.2 permet de conclure que $(A_3^N(\lfloor Nt \rfloor))_{t \ge 0}$ converge en loi vers le processus $(A_3(t))_{t \ge 0}$ quand $N$ tend vers l'infini. En d'autres termes, lorsque $N$ tend vers l'infini et que le temps est mesuré en unités de $N$ générations, le mécanisme d'apparition de l'ACPR d'un groupe composé de trois individus distincts pour le processus limite est le suivant :

- après un temps exponentiel $T_3 \sim \text{Exp}(3)$, un ancêtre commun à 2 individus apparaît,

- l'ancêtre commun apparaît alors après un temps $T_2 \sim \text{Exp}(1)$,

- les variables aléatoires $T_2$ et $T_3$ sont indépendantes.

En particulier, à la limite, la probabilité que le premier ancêtre commun soit commun aux trois individus est nulle. On généralise le résultat comme suit.

**Théorème 13.3** (Du discret au continu). *Pour tout $n \ge 2$, la convergence de processus suivante a lieu, au sens des lois marginales,*

$$ (A_n^N(\lfloor Nt \rfloor))_{t \ge 0} \xrightarrow[N \to \infty]{\text{loi}} (A_n(t))_{t \ge 0}. $$
Démonstration. On esquisse seulement les grandes lignes. Considérons dans un premier temps le comportement de la matrice de transition de $G_N$ lorsque $N$ est grand. Puisque $S_{k,k-1} = \binom{k}{2}$, pour tout $2 \le k \le n$,

$$g_{k,k-1}^{(N)} = S_{k,k-1} \frac{N(N-1)\cdots(N-k+1)}{N^k} = \binom{k}{2} \frac{1}{N} + \mathcal{O}(N^{-2})$$

tandis que pour $j < k - 1$, on a

$$g_{k,j}^{(N)} = S_{k,j} \frac{N(N-1)\cdots(N-j+1)}{N^k} = \mathcal{O}(N^{-2}).$$

Enfin, on a

$$g_{k,k}^{(N)} = N^{-k}N(N-1)\cdots(N-k+1) = 1 - \binom{k}{2} \frac{1}{N} + \mathcal{O}(N^{-2}).$$

Ainsi $G_N = I + N^{-1}Q + \mathcal{O}(N^{-2})$. La suite de la démonstration, omise ici, consiste à identifier la limite de $(G_N)^{[Nt]}$. $\square$

## 13.2 Longueur de l'arbre généalogique

Le processus à temps continu $(A_n(t))_{t \ge 0}$ d'espace d'états $\{1, \dots, n\}$ est un processus de mort pur, c'est-à-dire que ses trajectoires continues par morceaux sont décroissantes. Il est issu de $A_n(0) = n$ et décroît uniquement avec des sauts d'amplitude $-1$. Pour tout $k \ge 2$, sachant que le processus est dans l'état $k$, le temps d'attente $T_k$ avant de passer dans l'état $k-1$ (seule transition possible) suit une loi $\text{Exp}((\frac{k}{2}))$. De plus, les variables aléatoires $(T_k)_{2 \le k \le n}$ sont indépendantes. Le temps d'apparition de l'ACPR $W_n$ s'écrit donc

$$W_n = T_n + \dots + T_2,$$

où les $(T_k)_{2 \le k \le n}$ sont les temps de séjours dans les états $2, \dots, n$. En particulier, on a donc

$$E(W_n) = 2 \sum_{k=2}^{n} \frac{1}{k(k-1)} = 2 \sum_{k=2}^{n} \left( \frac{1}{k-1} - \frac{1}{k} \right) = 2 \left( 1 - \frac{1}{n} \right) \to 2$$

et les $T_2, \dots, T_n$ étant indépendantes, on a de plus

$$\begin{align*} \operatorname{Var}(W_n) &= \sum_{k=2}^{n} \operatorname{Var}(T_k) \\ &= \sum_{k=2}^{n} \frac{4}{k^2(k-1)^2} \end{align*}$$