Esercizio - Catena di Markov a tre stati, potenze di P e distribuzione stazionaria

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. Un server è in uno di tre stati a ogni secondo: $1=$ libero, $2=$ occupato, $3=$ guasto. Le transizioni, indipendenti dal passato, hanno matrice $$P=\begin{pmatrix}0{,}7&0{,}2&0{,}1\\0{,}4&0{,}5&0{,}1\\0{,}5&0&0{,}5\end{pmatrix}.$$

1. Verificare che $P$ è stocastica e calcolare $P^2$. Il server parte libero: probabilità di essere guasto dopo $3$ secondi.
2. Classificare la catena e calcolare la distribuzione stazionaria $\pi$.
3. Tempo medio di ritorno in ciascuno stato.
4. Tempo medio per arrivare al guasto partendo da libero e da occupato.
5. Quanti secondi servono perché $P^n$ approssimi $\pi$ a meno di $10^{-3}$ per entrata? Qual è la frazione di tempo guasto su $1000$ secondi?

Soluzione

1. Le righe sommano a $1$: $0{,}7+0{,}2+0{,}1=0{,}4+0{,}5+0{,}1=0{,}5+0+0{,}5=1$. Il prodotto riga per colonna dà, per esempio, $(P^2)_{11}=0{,}7\cdot0{,}7+0{,}2\cdot0{,}4+0{,}1\cdot0{,}5=0{,}62$: $$P^2=\begin{pmatrix}0{,}62&0{,}24&0{,}14\\0{,}53&0{,}33&0{,}14\\0{,}60&0{,}10&0{,}30\end{pmatrix}.$$ La probabilità di guasto al terzo secondo è $(P^3)_{13}$, cioè (riga 1 di $P^2$) per (colonna 3 di $P$): $0{,}62\cdot0{,}1+0{,}24\cdot0{,}1+0{,}14\cdot0{,}5=0{,}062+0{,}024+0{,}07=\boxed{0{,}156}$ $(=39/250)$.

2. Da $1$ si va in $2$ e in $3$; da $2$ in $1$; da $3$ in $1$: ogni stato raggiunge ogni altro (per esempio $2\to1\to3$ e $3\to1\to2$): irriducibile. Ha cappi ($P_{11}>0$), quindi è aperiodica. Per il teorema fondamentale esiste un'unica $\pi$, ed è il limite di $P^n$. Le equazioni $\pi P=\pi$ per le colonne $2$ e $3$:

• colonna 2: $\pi_2=0{,}2\pi_1+0{,}5\pi_2\Rightarrow0{,}5\pi_2=0{,}2\pi_1\Rightarrow\pi_2=0{,}4\pi_1$;
• colonna 3: $\pi_3=0{,}1\pi_1+0{,}1\pi_2+0{,}5\pi_3\Rightarrow0{,}5\pi_3=0{,}1\pi_1+0{,}04\pi_1=0{,}14\pi_1\Rightarrow\pi_3=0{,}28\pi_1$.

La colonna 1 è ridondante. Normalizzando: $\pi_1(1+0{,}4+0{,}28)=1{,}68\,\pi_1=1$, quindi $$\boxed{\pi=\left(\tfrac{25}{42},\ \tfrac5{21},\ \tfrac16\right)\approx(0{,}5952;\ 0{,}2381;\ 0{,}1667).}$$ Controllo con la colonna 1: $0{,}7\cdot\tfrac{25}{42}+0{,}4\cdot\tfrac{10}{42}+0{,}5\cdot\tfrac7{42}=\tfrac{17{,}5+4+3{,}5}{42}=\tfrac{25}{42}$.

3. Dal teorema $h_{ii}=1/\pi_i$: $h_{11}=\tfrac{42}{25}=1{,}68$, $h_{22}=\tfrac{21}5=4{,}2$, $h_{33}=6$ secondi. Partendo da guasto, il server torna guasto in media dopo $6$ secondi.

4. Sia $g_i$ il tempo medio per raggiungere lo stato $3$ da $i$. Condizionando al primo passo ($g_3=0$): $$g_1=1+0{,}7g_1+0{,}2g_2,\qquad g_2=1+0{,}4g_1+0{,}5g_2.$$ Dalla seconda $0{,}5g_2=1+0{,}4g_1$, quindi $g_2=2+0{,}8g_1$; nella prima: $0{,}3g_1=1+0{,}2(2+0{,}8g_1)=1{,}4+0{,}16g_1$, cioè $0{,}14g_1=1{,}4$, $g_1=10$ e $g_2=2+8=10$. In entrambi i casi $\boxed{10\text{ secondi}}$. Coerenza con il punto 3: da $3$ si resta (probabilità $0{,}5$, un passo) oppure si va in $1$ (probabilità $0{,}5$) e poi servono in media $g_1$ passi per tornare in $3$, quindi $h_{33}=1+0{,}5\,g_1=1+5=6$, come sopra.

5. Gli autovalori di $P$ sono $1$, $0{,}4$, $0{,}3$ (calcolati con numpy; la traccia $0{,}7+0{,}5+0{,}5=1{,}7$ è la loro somma). L'errore decade come $0{,}4^n$ a meno di una costante; le entrate della prima riga sono (potenze esatte):

Con $\pi_1=0{,}59524$ e $\pi_3=0{,}16667$, l'errore massimo su tutte le $9$ entrate di $P^n$ vale $2{,}4\cdot10^{-3}$ per $n=7$ e $9{,}99\cdot10^{-4}<10^{-3}$ per $n=8$ (calcolo numerico): servono $\boxed{n=8}$ secondi. La frazione di tempo guasto su $1000$ secondi è, per il teorema ergodico, $\approx\pi_3=\tfrac16$, cioè circa $167$ secondi.

Codice di verifica

import numpy as np
P = np.array([[.7, .2, .1], [.4, .5, .1], [.5, 0, .5]])
print(np.linalg.matrix_power(P, 2))
print("P^3[0,2] =", np.linalg.matrix_power(P, 3)[0, 2])
w, v = np.linalg.eig(P.T)
pi = np.real(v[:, np.argmax(np.real(w))]); pi /= pi.sum()
print("pi =", pi, " ritorni =", 1 / pi, " autovalori =", np.sort(np.real(w)))
# tempi di raggiungimento dello stato 3 (indice 2)
Q = P[:2, :2]
print("g =", np.linalg.solve(np.eye(2) - Q, np.ones(2)))
rng = np.random.default_rng(1); x = 0; n3 = 0
for _ in range(1000):
    x = rng.choice(3, p=P[x]); n3 += (x == 2)
print("secondi guasti su 1000:", n3)

Output (il numero di secondi guasti dipende dal seme):

[[0.62 0.24 0.14]
 [0.53 0.33 0.14]
 [0.6  0.1  0.3 ]]
P^3[0,2] = 0.156
pi = [0.5952381  0.23809524 0.16666667]  ritorni = [1.68 4.2  6.  ]  autovalori = [0.3 0.4 1. ]
g = [10. 10.]