Trouvez la loi stationnaire

Au chapitre 1, vous avez calculé la météo dans trois jours. Mais quelle est la proportion de jours de pluie sur une année entière ? Et cette proportion dépend-elle du temps qu’il faisait le 1er janvier ? L’intuition dit que non : au bout d’un moment, la chaîne « oublie » son point de départ.

Ce chapitre donne un sens précis à cette intuition. Vous allez calculer la loi stationnaire d’une chaîne, apprendre dans quels cas la chaîne converge vers elle, et à quelle vitesse. Nous reprenons la météo rennaise du chapitre 1, et nous finirons par l’algorithme qui a fait la fortune de Google.

La chaîne oublie son point de départ

Rappelons la matrice de la météo (soleil, nuages, pluie) :

$$P = \begin{pmatrix} 0{,}6 & 0{,}3 & 0{,}1 \\ 0{,}3 & 0{,}4 & 0{,}3 \\ 0{,}1 & 0{,}4 & 0{,}5 \end{pmatrix}$$

Calculons $p(n) = p(0)\,P^n$ pour deux points de départ opposés : un jour de soleil et un jour de pluie.

$n$	Départ au soleil	Départ sous la pluie
1	(0,600 ; 0,300 ; 0,100)	(0,100 ; 0,400 ; 0,500)
5	(0,357 ; 0,363 ; 0,280)	(0,335 ; 0,368 ; 0,297)
10	(0,346 ; 0,365 ; 0,288)	(0,346 ; 0,365 ; 0,289)
20	(0,346 ; 0,365 ; 0,288)	(0,346 ; 0,365 ; 0,288)

Au bout d’une dizaine de jours, les deux colonnes sont identiques à trois décimales : la loi de $X_n$ ne dépend plus de $X_0$. Elle s’est stabilisée sur un vecteur, qu’on appelle la loi stationnaire. Pourquoi « stationnaire » ? Parce que si $p(n)$ ne bouge plus, alors $p(n+1) = p(n)\,P = p(n)$.

Quel que soit le temps du premier jour, les deux courbes se confondent en une dizaine de jours sur la même ligne pointillée : la loi stationnaire Π.

Définir et calculer la loi stationnaire

Une loi stationnaire (ou invariante) est un vecteur ligne $\Pi = (\Pi_i)_{i \in E}$ qui vérifie

$$\Pi\, P = \Pi, \qquad \sum_{i \in E} \Pi_i = 1, \qquad \Pi_i \geq 0.$$

La première équation dit que si la chaîne suit la loi $\Pi$ à un instant, elle la suit encore à l’instant suivant, et donc pour toujours. La deuxième en fait une loi de probabilité. En langage d’algèbre linéaire, $\Pi$ est un vecteur propre à gauche de $P$ pour la valeur propre 1, normalisé pour sommer à 1.

Le système $\Pi P = \Pi$ contient une équation redondante (les lignes de $P$ sommant à 1, l’une se déduit des autres). On en supprime donc une, et on la remplace par la normalisation.

Exemple

Pour la météo, $\Pi P = \Pi$ donne, colonne par colonne :

• soleil : $\Pi_0 = 0{,}6\,\Pi_0 + 0{,}3\,\Pi_1 + 0{,}1\,\Pi_2$ ;
• pluie : $\Pi_2 = 0{,}1\,\Pi_0 + 0{,}3\,\Pi_1 + 0{,}5\,\Pi_2$.

La première s’écrit $0{,}4\,\Pi_0 = 0{,}3\,\Pi_1 + 0{,}1\,\Pi_2$, la seconde $0{,}5\,\Pi_2 = 0{,}1\,\Pi_0 + 0{,}3\,\Pi_1$. En les soustrayant, $0{,}5\,\Pi_0 = 0{,}6\,\Pi_2$, donc $\Pi_2 = \tfrac{5}{6}\Pi_0$. En reportant dans la première, $0{,}3\,\Pi_1 = 0{,}4\,\Pi_0 - \tfrac{0{,}5}{6}\Pi_0$, soit $\Pi_1 = \tfrac{19}{18}\Pi_0$. La normalisation donne enfin $\Pi_0\,(1 + \tfrac{19}{18} + \tfrac{15}{18}) = 1$, d’où

$$\Pi = \left(\tfrac{18}{52},\ \tfrac{19}{52},\ \tfrac{15}{52}\right) \approx (0{,}346 \,;\, 0{,}365 \,;\, 0{,}288).$$

Ce sont exactement les valeurs du tableau à $n = 20$. À long terme, il pleut environ 29 % des jours.

Pour une chaîne à deux états, le calcul se fait de tête. Reprenons la machine du chapitre 2, qui tombe en panne avec probabilité $a$ et est réparée avec probabilité $b$ :

$$P = \begin{pmatrix} 1-a & a \\ b & 1-b \end{pmatrix} \quad\Longrightarrow\quad \Pi = \left(\frac{b}{a+b},\ \frac{a}{a+b}\right).$$

En effet, la première équation donne $a\,\Pi_0 = b\,\Pi_1$ : le « flux » de la marche vers la panne égale le flux de la panne vers la marche. Avec $a = 0{,}1$ et $b = 0{,}4$, la machine fonctionne 80 % du temps.

Quand la loi stationnaire est-elle la limite ?

Une loi stationnaire existe-t-elle toujours ? Est-elle unique ? Et surtout, $p(n)$ converge-t-elle vers elle ? Le résultat central du chapitre répond à tout.

À retenir

Théorème ergodique. Si la chaîne est irréductible, récurrente positive (automatique si elle est finie) et apériodique, alors :

• il existe une unique loi stationnaire $\Pi$, solution de $\Pi P = \Pi$, $\sum_i \Pi_i = 1$ ;
• pour tous états $i$ et $j$, $\lim_{n \to \infty} p_{ij}(n) = \Pi_j$, et donc $\lim_{n \to \infty} p(n) = \Pi$ quelle que soit la loi initiale.

Chaque hypothèse compte. Voyons ce qui se passe quand l’une manque.

Sans apériodicité. Le robot du chapitre 2 fait la navette entre deux entrepôts : $P = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$. La loi $\Pi = (\tfrac12, \tfrac12)$ est bien stationnaire et unique. Mais partant de l’entrepôt 0, $p(n)$ vaut $(1, 0)$ aux instants pairs et $(0, 1)$ aux instants impairs : elle oscille et ne converge pas. En revanche, $\Pi$ garde un sens : c’est la proportion du temps passée dans chaque état sur une longue période (le robot est la moitié du temps dans chaque entrepôt).

Chaîne périodique : p₀(n) saute sans fin de 1 à 0 et ne converge pas, mais la proportion du temps passé en 0 converge bien vers Π₀ = 1/2.

Sans irréductibilité. Dans la chaîne du site du chapitre 2, deux classes sont fermées. La loi qui met toute la masse sur Départ est stationnaire, celle qui se répartit sur Abonné et Vidéo en proportions $(\tfrac{6}{11}, \tfrac{5}{11})$ aussi, et tout mélange des deux également. Il y a une infinité de lois stationnaires, et la limite dépend du point de départ. C’est l’objet du chapitre 4.

Sans récurrence positive (espace infini seulement). Dans une file d’attente où les clients arrivent plus vite qu’ils ne sont servis, la file grandit sans fin : aucune loi stationnaire n’existe. Nous le verrons au chapitre 8.

Piège

« Stationnaire » ne veut pas dire « limite ». Une chaîne périodique a une loi stationnaire vers laquelle elle ne converge pas. Avant d’écrire $\lim p(n) = \Pi$, vérifiez les trois hypothèses : irréductible, récurrente positive, apériodique.

Vitesse de convergence et temps de retour

À quelle vitesse la chaîne oublie-t-elle son départ ? Pour la machine, on peut calculer $p(n)$ exactement. Partant de l’état « marche » :

$$p_0(n) = \frac{b}{a+b} + \frac{a}{a+b}\,(1 - a - b)^n$$

Le premier terme est $\Pi_0$ ; le second est l’écart à la limite, qui décroît géométriquement avec la raison $1 - a - b$. Avec $a = 0{,}1$ et $b = 0{,}4$, on a $p_0(n) = 0{,}8 + 0{,}2 \times 0{,}5^n$ : l’écart est divisé par deux chaque heure (1 ; 0,9 ; 0,85 ; 0,825…).

Ce nombre $1 - a - b$ n’est pas un hasard : c’est la seconde valeur propre de $P$ (la première vaut toujours 1). En général, la vitesse de convergence est gouvernée par la plus grande valeur propre de $P$ en module après 1. Pour la météo, elle vaut environ 0,456 : l’écart à la loi stationnaire est divisé par un peu plus de 2 chaque jour, ce qui explique la stabilisation en une dizaine de jours. Une valeur proche de 1 signale une chaîne qui « se mélange » lentement.

L'écart à la limite est divisé par deux à chaque pas pour la machine, et par environ 2,2 pour la météo. En échelle logarithmique, une décroissance géométrique devient une droite, dont la pente est fixée par la seconde valeur propre.

La loi stationnaire a une autre interprétation, très parlante. Pour une chaîne irréductible et récurrente positive, la probabilité stationnaire d’un état est l’inverse de son temps moyen de retour :

$$\Pi_i = \frac{1}{\mu_i}$$

où $\mu_i$ est le nombre moyen de pas pour revenir en $i$ en partant de $i$. Plus on revient vite, plus on passe de temps dans l’état. Pour la machine, $\mu_0 = 1{,}25$ heure (calculé au chapitre 2), et on retrouve bien $\Pi_0 = 1/1{,}25 = 0{,}8$. Pour la météo, le temps moyen entre deux jours de pluie vaut $52/15 \approx 3{,}5$ jours.

PageRank : la loi stationnaire d’un surfeur aléatoire

Comment classer des pages web par importance ? L’idée de PageRank est de lâcher un surfeur aléatoire sur le web : sur chaque page, il clique au hasard sur l’un des liens sortants. Une page importante est une page où il passe souvent, c’est-à-dire une page de forte probabilité stationnaire.

Problème : le graphe du web n’est ni irréductible (des pages sans lien entrant, des culs-de-sac) ni forcément apériodique. On corrige avec la téléportation : à chaque pas, avec probabilité $1 - \alpha$, le surfeur saute sur une page tirée uniformément. La matrice devient

$$G = \alpha\, M + (1 - \alpha)\,\frac{1}{N}\,J$$

où $M$ est la matrice de la marche sur les liens, $N$ le nombre de pages, $J$ la matrice remplie de 1, et $\alpha$ le facteur d’amortissement, classiquement 0,85. Tous les coefficients de $G$ sont strictement positifs : la chaîne est irréductible et apériodique, et le théorème ergodique s’applique.

Exemple

Quatre pages : la page 1 pointe vers 2 et 3, la page 2 vers 3, la page 3 vers 1, la page 4 vers 3. Avec $\alpha = 0{,}85$, la loi stationnaire vaut environ $(0{,}373 \,;\, 0{,}196 \,;\, 0{,}394 \,;\, 0{,}038)$. La page 3, qui reçoit trois liens, arrive en tête ; la page 4, qui n’en reçoit aucun, n’a que la part due à la téléportation, $0{,}15/4 = 0{,}0375$. La page 1 est presque aussi bien classée que la page 3 alors qu’elle ne reçoit qu’un lien : il vient de la page la plus importante.

La page 3 reçoit trois liens et arrive en tête ; la page 1, qui n'en reçoit qu'un mais depuis la page 3, la suit de près. La page 4, sans lien entrant, n'a que la part de la téléportation (pointillés).

Pour calculer $\Pi$ sur des milliards de pages, on ne résout pas de système linéaire : on itère simplement $\pi \leftarrow \pi\, G$ à partir d’une loi quelconque. C’est la méthode de la puissance, qui converge d’autant plus vite que la seconde valeur propre est petite (pour PageRank, elle est au plus $\alpha$).

import numpy as np

P = np.array([[0.6, 0.3, 0.1],
              [0.3, 0.4, 0.3],
              [0.1, 0.4, 0.5]])
n = len(P)

# 1) Système linéaire : (P^T - I) Pi^T = 0, une équation remplacée par sum(Pi) = 1
A = P.T - np.eye(n)
A[-1, :] = 1
b = np.zeros(n); b[-1] = 1
pi = np.linalg.solve(A, b)
print(pi)                                   # [0.3462 0.3654 0.2885]

# 2) Méthode de la puissance
pi_it = np.ones(n) / n
for _ in range(100):
    pi_it = pi_it @ P
print(np.allclose(pi, pi_it))               # True

# 3) Valeurs propres : 1, puis la seconde qui fixe la vitesse de convergence
print(np.sort(np.abs(np.linalg.eigvals(P)))[::-1])   # [1. 0.456 0.044]

# 4) Proportion du temps sur une longue trajectoire simulée
rng = np.random.default_rng(1)
x, compte = 0, np.zeros(n)
for _ in range(100_000):
    x = rng.choice(n, p=P[x])
    compte[x] += 1
print(compte / compte.sum())                # proche de pi

Le quatrième calcul illustre une autre face du théorème ergodique : la proportion du temps passée dans chaque état le long d’une seule trajectoire converge elle aussi vers $\Pi$. C’est ce résultat qui fonde les méthodes MCMC en statistique bayésienne : on construit une chaîne dont la loi stationnaire est la loi qu’on veut échantillonner, on la laisse tourner un moment (la « chauffe », qui correspond au régime transitoire), puis on utilise ses états comme des tirages.

Une seule trajectoire suffit : au fil des jours, la fraction du temps passée dans chaque état se cale sur la loi stationnaire (pointillés). Les premières centaines de jours, encore marqués par le départ, jouent le rôle de la chauffe en MCMC.

En entretien

« La loi stationnaire est le vecteur propre à gauche de la matrice de transition pour la valeur propre 1, normalisé. Elle est la limite de la chaîne si celle-ci est irréductible, récurrente positive et apériodique ; PageRank ajoute la téléportation précisément pour garantir ces conditions. »

En résumé

• Une loi stationnaire vérifie $\Pi P = \Pi$, $\sum_i \Pi_i = 1$ : c’est un vecteur propre à gauche de $P$ pour la valeur propre 1.
• On la calcule en remplaçant une équation de $\Pi P = \Pi$ par la normalisation, ou par la méthode de la puissance.
• Théorème ergodique : irréductible + récurrente positive (automatique si finie) + apériodique implique une unique $\Pi$ et $p(n) \to \Pi$ quelle que soit la loi initiale.
• Périodique : $\Pi$ existe mais $p(n)$ oscille ; $\Pi$ reste la proportion du temps passée dans chaque état. Non irréductible : plusieurs lois stationnaires possibles.
• La convergence est géométrique, à la vitesse de la seconde plus grande valeur propre en module.
• $\Pi_i = 1/\mu_i$ : probabilité stationnaire et temps moyen de retour sont inverses. PageRank est la loi stationnaire d’une marche aléatoire avec téléportation.

Exercices

Exercice · Disponibilité d'une machine

facile

Une imprimante 3D tombe en panne en cours d’heure avec probabilité 0,2 et, quand elle est en panne, est réparée dans l’heure avec probabilité 0,3. Quelle proportion du temps fonctionne-t-elle à long terme ? En moyenne, combien d’heures s’écoulent entre deux heures de panne ?

Voir la solution

Avec $a = 0{,}2$ et $b = 0{,}3$, la chaîne est irréductible, finie et apériodique (boucles), donc ergodique. La loi stationnaire vaut

$$\Pi = \left(\frac{b}{a+b},\ \frac{a}{a+b}\right) = (0{,}6 \,;\, 0{,}4).$$

L’imprimante fonctionne 60 % du temps. Le temps moyen de retour dans l’état « panne » vaut $\mu_1 = 1/\Pi_1 = 2{,}5$ heures.

En code

Python

import numpy as np

a, b = 0.2, 0.3                          # panne, réparation (par heure)
P = np.array([[1 - a, a],
              [b, 1 - b]])

# Pi P = Pi : une équation remplacée par la normalisation
A = P.T - np.eye(2)
A[-1, :] = 1
pi = np.linalg.solve(A, [0, 1])
print(pi)                                # [0.6 0.4] : en marche 60 % du temps
print(1 / pi[1])                         # 2.5 heures entre deux heures de panne

R

a <- 0.2; b <- 0.3                      # panne, réparation (par heure)
P <- matrix(c(1 - a, a,
              b, 1 - b), nrow = 2, byrow = TRUE)

# Pi P = Pi : une équation remplacée par la normalisation
A <- t(P) - diag(2)
A[2, ] <- 1
pi_ <- solve(A, c(0, 1))
pi_                                     # 0.6 0.4 : en marche 60 % du temps
1 / pi_[2]                              # 2.5 heures entre deux heures de panne

Exercice · Une chaîne à trois états

moyen

Un livreur travaille dans trois quartiers alignés A, B, C. Depuis A, il reste en A ou passe en B avec probabilité 1/2 chacun ; depuis B, il va en A ou en C avec probabilité 1/4 chacun et reste en B avec probabilité 1/2 ; depuis C, il reste en C ou passe en B avec probabilité 1/2 chacun.

1. Justifiez que la chaîne converge vers une loi stationnaire unique.
2. Calculez-la, puis le temps moyen entre deux passages en A.

Voir la solution

1. Tous les quartiers communiquent (A ⇄ B ⇄ C) : la chaîne est irréductible. Elle est finie, donc récurrente positive, et chaque état a une boucle, donc elle est apériodique. Le théorème ergodique s’applique.
2. L’équation de la colonne A donne $\Pi_A = \tfrac12 \Pi_A + \tfrac14 \Pi_B$, soit $\Pi_B = 2\,\Pi_A$. Par symétrie entre A et C, $\Pi_C = \Pi_A$. La normalisation donne $\Pi_A (1 + 2 + 1) = 1$, donc

$$\Pi = \left(\tfrac14,\ \tfrac12,\ \tfrac14\right).$$

Le temps moyen de retour en A vaut $\mu_A = 1/\Pi_A = 4$ pas.

En code

Python

import numpy as np

P = np.array([[1/2, 1/2, 0],             # quartiers A, B, C
              [1/4, 1/2, 1/4],
              [0, 1/2, 1/2]])

# Vecteur propre à gauche pour la valeur propre 1, normalisé
valeurs, vecteurs = np.linalg.eig(P.T)
pi = np.real(vecteurs[:, np.argmin(abs(valeurs - 1))])
pi = pi / pi.sum()
print(pi)                                # [0.25 0.5  0.25]
print(1 / pi[0])                         # 4.0 pas entre deux passages en A
print(np.linalg.matrix_power(P, 50)[2])  # [0.25 0.5  0.25] : partant de C, on converge

R

P <- matrix(c(1/2, 1/2, 0,              # quartiers A, B, C
              1/4, 1/2, 1/4,
              0,   1/2, 1/2), nrow = 3, byrow = TRUE)

# Vecteur propre à gauche pour la valeur propre 1, normalisé
e <- eigen(t(P))
pi_ <- Re(e$vectors[, which.min(abs(e$values - 1))])
pi_ <- pi_ / sum(pi_)
pi_                                     # 0.25 0.50 0.25
1 / pi_[1]                              # 4 pas entre deux passages en A

Exercice · Stationnaire mais pas limite

difficile

Une balle rebondit entre trois positions alignées : depuis une extrémité (0 ou 2), elle va toujours au centre (1) ; depuis le centre, elle va à chaque extrémité avec probabilité 1/2.

1. Quelle est la période de la chaîne ? Calculez sa loi stationnaire.
2. La balle part de la position 0. Calculez $p(1)$, $p(2)$, $p(3)$. La suite $p(n)$ converge-t-elle ?
3. Montrez que la moyenne de $p(n)$ sur deux pas consécutifs est égale à la loi stationnaire. Interprétez.

Voir la solution

1. Tout chemin alterne entre le centre et une extrémité, donc un retour n’est possible qu’en un nombre pair de pas : la période vaut 2. La colonne 0 de $\Pi P = \Pi$ donne $\Pi_0 = \tfrac12 \Pi_1$, de même $\Pi_2 = \tfrac12 \Pi_1$, et la normalisation conduit à $\Pi = (\tfrac14, \tfrac12, \tfrac14)$.
2. $p(0) = (1, 0, 0)$, $p(1) = (0, 1, 0)$, $p(2) = (\tfrac12, 0, \tfrac12)$, $p(3) = (0, 1, 0)$, puis la suite alterne entre $(\tfrac12, 0, \tfrac12)$ et $(0, 1, 0)$. Elle ne converge pas : la chaîne est irréductible et récurrente positive, mais périodique.
3. $\tfrac12\big[(\tfrac12, 0, \tfrac12) + (0, 1, 0)\big] = (\tfrac14, \tfrac12, \tfrac14) = \Pi$. La loi stationnaire ne décrit pas la position à un instant donné, mais la proportion du temps passée à chaque position : la balle est au centre un pas sur deux, et à chaque extrémité un pas sur quatre.

En code

Python

import numpy as np

P = np.array([[0, 1, 0],                 # extrémité 0 -> centre
              [0.5, 0, 0.5],             # centre -> une extrémité
              [0, 1, 0]])
p = np.array([1.0, 0, 0])                # départ en 0
lois = [p]
for n in range(5):
    p = p @ P
    lois.append(p)
for n, loi in enumerate(lois):
    print(n, loi)
# 0 [1. 0. 0.]
# 1 [0. 1. 0.]
# 2 [0.5 0.  0.5]
# 3 [0. 1. 0.]      et ainsi de suite : la suite oscille
print((lois[4] + lois[5]) / 2)           # [0.25 0.5  0.25] = Pi

R

P <- matrix(c(0,   1, 0,                # extrémité 0 -> centre
              0.5, 0, 0.5,              # centre -> une extrémité
              0,   1, 0), nrow = 3, byrow = TRUE)
p <- c(1, 0, 0)                         # départ en 0
lois <- list(p)
for (n in 1:5) {
  p <- as.vector(p %*% P)
  lois[[n + 1]] <- p
}
do.call(rbind, lois)                    # lignes : n = 0, ..., 5
## 1.0 0 0.0 / 0 1 0 / 0.5 0 0.5 / 0 1 0 / 0.5 0 0.5 / 0 1 0 : la suite oscille
(lois[[5]] + lois[[6]]) / 2             # 0.25 0.50 0.25 = Pi