Comment simuler une chaîne de… : En itérant sur numpy.random.choice(range(… — ECG2 appliquées | MetMat

Comment simuler une chaîne de Markov et observer son état stable en Python ?

En itérant sur numpy.random.choice(range(r), p=M[etat_courant]) et en observant la convergence vers l'état stable

Approfondissement — Simuler la trajectoire aléatoire d'une chaîne de Markov sur un graphe fini et estimer empiriquement sa loi stationnaire.

Maîtriser

Reconnaître

Résoudre

Hors programme — Cette méthode va au-delà du B.O. officiel. Proposée pour aller plus loin.

Que ferais-tu face à un exercice de ce type ?

Météo simple : trois états $\{S, N, P\}$ (Soleil, Nuageux, Pluie) avec matrice de transition $M = \begin{pmatrix} 0{,}7 & 0{,}2 & 0{,}1 \\ 0{,}3 & 0{,}4 & 0{,}3 \\ 0{,}2 & 0{,}3 & 0{,}5 \end{pmatrix}$ . Estimer la loi stationnaire par simulation sur $N = 10\,000$ jours.

L'objectif

Approfondissement — Simuler la trajectoire aléatoire d'une chaîne de Markov sur un graphe fini et estimer empiriquement sa loi stationnaire.

Le principe

Si $M$ est la matrice de transition (ligne $i$ = loi conditionnelle sachant l'état $i$ ), alors numpy.random.choice(range(r), p=M[i]) tire l'état suivant avec la bonne probabilité ; pour une chaîne irréductible apériodique, la fréquence empirique des états converge vers l'unique état stable $\pi$ vérifiant $\pi M = \pi$ .

La méthode

1
Je code la matrice de transition $M$ en tableau numpy de taille $r \times r$ , en vérifiant que chaque ligne somme à $1$ .
2
Je fixe l'état initial etat = i0 et le nombre N d'itérations, puis j'initialise un tableau trajectoire de taille N+1.
3
Pour n de 0 à N-1, je tire l'état suivant par etat = numpy.random.choice(range(r), p=M[etat]), que je stocke dans la trajectoire.
4
J'estime la loi stationnaire par fréquences : freq[k] = (trajectoire == k).mean(), et je compare au vecteur propre à gauche $\pi$ associé à la valeur propre $1$ de $M$ .

Pour comprendre, fais des exercices !

Exercice d'exemple

Étape 1 —

Je code la matrice de transition

M

en tableau numpy de taille

r \times r

, en vérifiant que chaque ligne somme à

1

Je code la matrice :

import numpy as np
M = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])
assert np.allclose(M.sum(axis=1), 1)

Étape 2 —

Je fixe l'état initial etat = i0 et le nombre N d'itérations, puis j'initialise un tableau trajectoire de taille N+1.

Je fixe l'état initial et la longueur :

N = 10_000
etat = 0  # on part de Soleil
traj = np.zeros(N+1, dtype=int)
traj[0] = etat

Étape 3 —

Pour n de 0 à N-1, je tire l'état suivant par etat = numpy.random.choice(range(r), p=M[etat]), que je stocke dans la trajectoire.

Je simule la chaîne :

for n in range(N):
    etat = np.random.choice(3, p=M[etat])
    traj[n+1] = etat

Étape 4 —

J'estime la loi stationnaire par fréquences : freq[k] = (trajectoire == k).mean(), et je compare au vecteur propre à gauche

\pi

associé à la valeur propre

1

M

Je calcule les fréquences empiriques :

freq = np.array([(traj == k).mean() for k in range(3)])
print(freq)  # ex: [0.45, 0.27, 0.28]

Les fréquences sont proches de l'état stable théorique $\pi \approx (0{,}45\,;\, 0{,}28\,;\, 0{,}27)$ .

L'état stable empirique vaut environ $(0{,}45\,;\, 0{,}28\,;\, 0{,}27)$ , en accord avec le calcul théorique de $\pi$ tel que $\pi M = \pi$ .

Les applications suivantes et la correction guidée sont réservées aux membres Premium

La méthode et sa première application corrigée restent en accès libre. Le Premium débloque les applications suivantes, l'aide IA et le suivi de ta maîtrise.