Esercizio - Dimostrazioni sulle catene di Markov - potenze, aperiodicità e matrici doppiamente stocastiche

Teoria: Catene di Markov - matrice di transizione, classificazione degli stati e distribuzione stazionariaUna catena di Markov a tempo discreto e stati finiti è descritta da una matrice stocastica P: P(X_{n+1}=j | X_n=i)=P_ij, la legge al tempo n è π_0 P^n. Gli stati si classificano in classi comunicanti (irriducibilità), periodo (aperiodicità) e ricorrenza/transitorietà. Una distribuzione stazionaria soddisfa πP=π; per catene finite irriducibili esiste ed è unica, vale π_i=1/h_ii (h_ii = tempo medio di primo ritorno), e se la catena è anche aperiodica P^n converge a una matrice con tutte le righe uguali a π. Le catene reversibili (π_i P_ij = π_j P_ji) hanno stazionaria immediata.Catene di Markov - matrice di transizione, classificazione degli stati e distribuzione stazionaria →.

Testo (domande da orale, con verifica numerica).

1. Per $P=\begin{pmatrix}1-a&a\\b&1-b\end{pmatrix}$ con $a+b>0$ dimostrare per induzione la formula chiusa di $P^n$. Applicarla a $a=0{,}1$, $b=0{,}5$ per calcolare $P^5$ e il limite.
2. Dimostrare che una catena irriducibile con almeno uno stato $i$ tale che $P_{ii}>0$ è aperiodica.
3. Una matrice è doppiamente stocastica se anche le colonne sommano a $1$. Dimostrare che la distribuzione uniforme è stazionaria. Applicazione: $P=\begin{pmatrix}0{,}2&0{,}5&0{,}3\\0{,}5&0{,}1&0{,}4\\0{,}3&0{,}4&0{,}3\end{pmatrix}$, probabilità di stare nello stato $1$ a lungo termine e tempo medio di ritorno in ciascuno stato.
4. Dimostrare che se $\pi$ è stazionaria e la catena parte da $\pi$, allora $P(X_0=i,X_1=j)=\pi_iP_{ij}$ e $X_1\sim\pi$; dedurre che nelle catene reversibili la legge del processo invertito nel tempo coincide con quella originale (verifica con la coppia $(X_0,X_1)$ nel modello $Q$ della teoria).

Soluzione

1. Poniamo $\Pi=\dfrac1{a+b}\begin{pmatrix}b&a\\b&a\end{pmatrix}$, $D=\dfrac1{a+b}\begin{pmatrix}a&-a\\-b&b\end{pmatrix}$ e $\lambda=1-a-b$. Tesi: $P^n=\Pi+\lambda^nD$.

• Base $n=1$. $\Pi+\lambda D=\dfrac1{a+b}\begin{pmatrix}b+\lambda a&a-\lambda a\\b-\lambda b&a+\lambda b\end{pmatrix}$. Si ha $b+(1-a-b)a=b+a-a^2-ab=(a+b)(1-a)$ e $a-\lambda a=a(a+b)$, quindi la prima riga è $(1-a,a)$; analogamente la seconda è $(b,1-b)$: è $P$.
• Due identità. (i) $\Pi P=\Pi$: la riga $(b,a)$ moltiplicata per $P$ dà $(b(1-a)+ab,\ ba+a(1-b))=(b,a)$ (è $\pi P=\pi$ con $\pi\propto(b,a)$). (ii) $DP=\lambda D$: $(a,-a)P=(a(1-a)-ab,\ a^2-a(1-b))=(a\lambda,-a\lambda)$, e analogamente la seconda riga.
• Passo. $P^{n+1}=P^nP=\Pi P+\lambda^nDP=\Pi+\lambda^{n+1}D$. $\blacksquare$

Per $a,b\in[0,1]$ con $a+b>0$ si ha $\lambda\in[-1,1)$, e $|\lambda|<1$ salvo $a=b=1$: allora $P^n\to\Pi$. Con $a=0{,}1$, $b=0{,}5$, $\lambda=0{,}4$: $\Pi$ ha righe $(\tfrac56,\tfrac16)$ e per $n=5$ si ha $\lambda^5=0{,}01024$, quindi $$P^5_{11}=\tfrac56+\tfrac{0{,}1}{0{,}6}\cdot0{,}01024=0{,}83333+0{,}0017067=0{,}83504,\qquad P^5_{12}=0{,}16496.$$ Si noti che per $\lambda=-1$ ($a=b=1$) non c'è convergenza: è il ciclo di periodo $2$.

2. Sia $P_{ii}>0$ e la catena irriducibile. Per $i$: $(P^1)_{ii}>0$, quindi $\gcd\{n:(P^n)_{ii}>0\}$ divide $1$, cioè $d(i)=1$. Il periodo è costante sulla classe (se $j\leftrightarrow i$ con cammini di lunghezze $r$ da $i$ a $j$ e $s$ da $j$ a $i$, ogni ritorno in $i$ di lunghezza $m$ dà ritorni in $j$ di lunghezze $r+s$ e $r+m+s$, quindi $d(j)$ divide la loro differenza $m$, e di conseguenza $d(j)\mid d(i)$; per simmetria sono uguali). L'unica classe è tutta la catena: aperiodica. $\blacksquare$

3. Sia $u=(\tfrac1k,\dots,\tfrac1k)$. La componente $j$ di $uP$ è $\tfrac1k\sum_iP_{ij}$, che vale $\tfrac1k$ perché la colonna $j$ somma a $1$. Quindi $uP=u$. $\blacksquare$

Applicazione: le colonne sommano a $0{,}2+0{,}5+0{,}3=1$, $0{,}5+0{,}1+0{,}4=1$, $0{,}3+0{,}4+0{,}3=1$ ✓. Tutte le entrate sono positive, quindi la catena è irriducibile e aperiodica: $\pi=(\tfrac13,\tfrac13,\tfrac13)$ è l'unica stazionaria, la probabilità di stare in $1$ a lungo termine è $\tfrac13$ e i tempi medi di ritorno sono $h_{ii}=1/\pi_i=3$ per ogni stato. Numericamente, la prima riga di $P^{10}$ è $(0{,}33335;\ 0{,}33331;\ 0{,}33334)$, già a meno di $5\cdot10^{-5}$ da $\tfrac13$.

4. $P(X_0=i,X_1=j)=P(X_0=i)P(X_1=j\mid X_0=i)=\pi_iP_{ij}$. Sommando su $i$ si ha $P(X_1=j)=\sum_i\pi_iP_{ij}=\pi_j$. Se vale anche il bilancio dettagliato, $P(X_0=i,X_1=j)=\pi_iP_{ij}=\pi_jP_{ji}=P(X_0=j,X_1=i)$: la coppia $(X_0,X_1)$ ha la stessa legge della coppia scambiata, cioè guardare il processo all'indietro non cambia le statistiche a due tempi. Per la catena $Q=\begin{pmatrix}\tfrac12&\tfrac12&0\\\tfrac14&\tfrac12&\tfrac14\\0&\tfrac12&\tfrac12\end{pmatrix}$ con $\pi=(\tfrac14,\tfrac12,\tfrac14)$ la tabella $\pi_iQ_{ij}$ è $$\begin{pmatrix}\tfrac18&\tfrac18&0\\\tfrac18&\tfrac14&\tfrac18\\0&\tfrac18&\tfrac18\end{pmatrix},$$ simmetrica (e con somma $1$): reversibile ✓.

Codice di verifica

import numpy as np
a, b = 0.1, 0.5
P = np.array([[1 - a, a], [b, 1 - b]])
lam = 1 - a - b
Pi = np.array([[b, a], [b, a]]) / (a + b); D = np.array([[a, -a], [-b, b]]) / (a + b)
for n in range(1, 8):
    assert np.allclose(np.linalg.matrix_power(P, n), Pi + lam**n * D)
print("P^5 =", np.linalg.matrix_power(P, 5)[0])
R = np.array([[.2, .5, .3], [.5, .1, .4], [.3, .4, .3]])
print("colonne:", R.sum(axis=0), " R^10 riga 1:", np.linalg.matrix_power(R, 10)[0])
Q = np.array([[.5, .5, 0], [.25, .5, .25], [0, .5, .5]]); pi = np.array([.25, .5, .25])
J = pi[:, None] * Q
print(J, " simmetrica:", np.allclose(J, J.T))

Output:

P^5 = [0.83504 0.16496]
colonne: [1. 1. 1.]  R^10 riga 1: [0.33335081 0.33330947 0.33333973]
[[0.125 0.125 0.   ]
 [0.125 0.25  0.125]
 [0.    0.125 0.125]]  simmetrica: True