import numpy as np
import matplotlib.pyplot as plt
import random
On part d'une matrice de transition $$ P= \begin{bmatrix} 0.2 & 0.7 & 0.1 \\ 0.9 & 0 & 0.1 \\ 0.2 & 0.8 & 0 \end{bmatrix} $$ et d'état initial $(1,0,0).$
# Matrice de transition
P = np.array([[0.2, 0.7, 0.1],
[0.9, 0. , 0.1],
[0.2, 0.8, 0.]])
# Itération de la chaîne
n=100 # nombre de pas
state = [1.0, 0.0, 0.0] # Etat initial
stateHist = [state] # Historique de tous les états
for x in range(n):
state = np.dot(state,P) #! Numpy convertit une liste en array pour np.dot
stateHist.append(state)
plt.plot(stateHist)
plt.show()
print(state)
Une façon de trouver la mesure invariante est de résoudre le système linéaire :
$(\,{}^tP-I)\pi=0$ sous la contrainte linéaire $\pi\cdot \boldsymbol 1=1$, ce que l'on peut écrire $A\pi=b$ avec $$ A:= \begin{bmatrix} P^t-I_3\\ \boldsymbol 1^t \end{bmatrix},\qquad b:=e_4. $$ Comme expliqué en cours, pour trouver $\pi$ il est equivalent de résoudre $A^t A \pi=A^tb$.
NB : On pourrait aussi minimiser $\|Ax-b\|^2_2$ (moindres carrés).
A = np.append(P.T-np.identity(3),[[1,1,1]],axis=0)
b = np.array([0,0,0,1])
#np.linalg.solve(np.transpose(A).dot(A),np.transpose(A).dot(b))
np.linalg.solve(np.dot(A.T,A),np.dot(A.T,b))
Enfin, si la chaîne satisfait la condition de Doeblin (il existe $n_0\geq 1$, $\beta>0$ et une mesure de probabilité $\nu$ sur $E$ tels que, pour tout $x,y\in E$, on a $(P^{n_0})_{xy}\geq \beta\nu_y$), alors pour tout $n\geq 1$, $$ \sum_{x\in E}|(P^n)_{xy}-\pi_x|\leq 2(1-\beta)^{\lfloor n/n_0\rfloor}. $$ On peut donc quantifier la vitesse de convergence de $P^n_{xy}$ vers $\pi_x$ en remarquant ici que la condition est satisfaite pour $n_0=1$, $\beta=0.2$ et $\nu=(1,0,0)$. Plus généralement, il suffit que $P^{n_0}$ ait une colonne strictement positive (et $\beta$ est l'inf sur cette colonne) pour que la condition de Doeblin soit satisfaite.
Simuler une chaine de Markov¶
Cela revient à tirer une multinomiale à chaque étape.
def Markov(n,mu,transition):
P = transition
state = random.choices(range(3),mu)
trajectory = state
for i in range(n):
p = P[state[0]]
state = random.choices(range(3),p) #! random.choices renvoie une liste
trajectory.append(state[0])
return trajectory
plt.plot(Markov(100,[1,0,0],P),
color='blue', marker='o', linestyle='dashed',
linewidth=1, markersize=4) # rajouter labels axes
plt.show()
