Esercizio - Modello di Ehrenfest con quattro palline e bilancio dettagliato

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 →; legge binomiale in Bernoulli, binomiale e uniforme discretaBernoulli Be(p): un solo tentativo, vale 1 con probabilità p; binomiale Bin(n, p): numero di successi in n prove indipendenti, P(X = k) = C(n,k) p^k (1−p)^(n−k), media np e varianza np(1−p); uniforme discreta: n valori equiprobabili.Bernoulli, binomiale e uniforme discreta →.

Testo. Quattro palline numerate stanno in due urne $A$ e $B$. A ogni passo si sceglie una pallina a caso (probabilità $\tfrac14$ ciascuna) e la si sposta nell'altra urna. Sia $X_n$ il numero di palline in $A$.

1. Scrivere la matrice di transizione sugli stati $0,\dots,4$.
2. Dimostrare con il bilancio dettagliato che $\pi=\text{Bin}(4,\tfrac12)$ è stazionaria.
3. Tempi medi di ritorno in $0$, $1$, $2$.
4. Periodo della catena. Calcolare la legge di $X_n$ partendo da $X_0=0$ per $n=1,\dots,6$ e per $n=20,21$.
5. Tempo medio per passare da $0$ a $2$ e da $0$ a $4$.
6. Con la catena pigra $\tfrac12(I+P)$ cosa cambia?

Soluzione

1. Se in $A$ ci sono $i$ palline, con probabilità $\tfrac i4$ si sceglie una pallina di $A$ (e $X$ scende a $i-1$) e con probabilità $\tfrac{4-i}4$ una di $B$ (e $X$ sale a $i+1$): $$P=\begin{pmatrix}0&1&0&0&0\\\tfrac14&0&\tfrac34&0&0\\0&\tfrac12&0&\tfrac12&0\\0&0&\tfrac34&0&\tfrac14\\0&0&0&1&0\end{pmatrix}.$$

2. Il bilancio dettagliato richiede $\pi_iP_{i,i+1}=\pi_{i+1}P_{i+1,i}$ per $i=0,\dots,3$. Con $\pi_i=\binom4i/16$: $$\frac{\binom4i}{16}\cdot\frac{4-i}4\quad\text{e}\quad\frac{\binom4{i+1}}{16}\cdot\frac{i+1}4$$ sono uguali perché $\binom4i(4-i)=\binom4{i+1}(i+1)$ (entrambi valgono $\tfrac{4!}{i!\,(3-i)!}$). Per $i=0$: $\tfrac1{16}\cdot1=\tfrac4{16}\cdot\tfrac14$ ✓. Dunque $$\pi=\left(\tfrac1{16},\tfrac4{16},\tfrac6{16},\tfrac4{16},\tfrac1{16}\right)=(0{,}0625;\ 0{,}25;\ 0{,}375;\ 0{,}25;\ 0{,}0625).$$ È l'interpretazione "ogni pallina è in $A$ con probabilità $\tfrac12$, indipendentemente": equilibrio che risulta naturale. La catena è irriducibile, quindi $\pi$ è l'unica stazionaria.

3. $h_{ii}=1/\pi_i$: $h_{00}=16$, $h_{11}=4$, $h_{22}=\tfrac{16}6=\tfrac83\approx2{,}67$ (e per simmetria $h_{33}=4$, $h_{44}=16$). Per tornare a "tutte le palline in $B$" ci vogliono in media $16$ passi, mentre lo stato $2$, tipico, si rivisita ogni $2{,}67$ passi.

4. A ogni passo $X$ cambia di $\pm1$, quindi la parità di $X_n$ alterna: da $0$ si torna a $0$ solo a tempi pari, il periodo è $2$ (grafo bipartito pari/dispari). La legge di $X_n$ da $X_0=0$ (calcolata con $\pi^{(n)}=\pi^{(n-1)}P$):

(Per esempio $n=2$: da $0$ si va in $1$ e da $1$ si va in $0$ con $\tfrac14$ o in $2$ con $\tfrac34$.) Per $n$ pari la legge tende a $(\tfrac18,0,\tfrac34,0,\tfrac18)$, per $n$ dispari a $(0,\tfrac12,0,\tfrac12,0)$: nessun limite, anche se $\pi$ esiste. Gli autovalori sono $\pm1,\pm\tfrac12,0$: l'autovalore $-1$ segnala il periodo $2$. La media ergodica $\tfrac12\bigl(\pi^{(20)}+\pi^{(21)}\bigr)=(\tfrac1{16},\tfrac14,\tfrac38,\tfrac14,\tfrac1{16})$ è proprio $\pi$.

5. Sia $e_i$ il tempo medio per passare da $i$ a $i+1$. Da $0$ si sale per forza: $e_0=1$. Da $i\ge1$ si sale con probabilità $\tfrac{4-i}4$ in un passo; con probabilità $\tfrac i4$ si scende a $i-1$ e poi bisogna tornare in $i$ (media $e_{i-1}$) e riprovare (media $e_i$): $$e_i=1+\tfrac i4\,(e_{i-1}+e_i)\ \Longrightarrow\ e_i\Bigl(1-\tfrac i4\Bigr)=1+\tfrac i4e_{i-1}.$$

• $e_1$: $\tfrac34e_1=1+\tfrac14\cdot1=\tfrac54$, $e_1=\tfrac53$.
• $e_2$: $\tfrac12e_2=1+\tfrac12\cdot\tfrac53=\tfrac{11}6$, $e_2=\tfrac{11}3$.
• $e_3$: $\tfrac14e_3=1+\tfrac34\cdot\tfrac{11}3=\tfrac{15}4$, $e_3=15$.

Quindi da $0$ a $2$: $e_0+e_1=1+\tfrac53=\boxed{\tfrac83\approx2{,}67}$; da $0$ a $4$: $1+\tfrac53+\tfrac{11}3+15=\boxed{\tfrac{64}3\approx21{,}33}$ passi. Salire fino alla configurazione estrema costa molto più che tornare allo stato tipico: lo stato $4$ ha probabilità stazionaria piccola ($\tfrac1{16}$) e la catena è "attratta" verso il centro (deriva verso $2$). Controllo: $h_{00}=1+e_{\,\text{da }1\text{ a }0}$, e per simmetria il tempo da $1$ a $0$ è $15$, quindi $h_{00}=1+15=16$ ✓.

6. La catena $\tilde P=\tfrac12(I+P)$ (a ogni passo, con probabilità $\tfrac12$ non si fa nulla) ha la stessa $\pi$ (da $\pi\tilde P=\tfrac12(\pi+\pi P)=\pi$) ed è aperiodica (cappi). Gli autovalori diventano $\tfrac{1+\lambda}2$: $1,\tfrac34,\tfrac12,\tfrac14,0$. Ora la legge converge: da $0$, a $n=20$ vale $(0{,}0633;\ 0{,}2516;\ 0{,}375;\ 0{,}2484;\ 0{,}0617)$, vicinissima a $\pi$; il passo di contrazione è $\tfrac34$ per passo.

Codice di verifica

import numpy as np
from math import comb
N = 4
P = np.zeros((N + 1, N + 1))
for i in range(N + 1):
    if i > 0: P[i, i - 1] = i / N
    if i < N: P[i, i + 1] = (N - i) / N
pi = np.array([comb(N, k) for k in range(N + 1)]) / 2**N
print("pi P == pi:", np.allclose(pi @ P, pi), " ritorni:", 1 / pi)
x = np.eye(N + 1)[0]
for n in range(1, 22):
    x = x @ P
    if n <= 6 or n >= 20: print(n, np.round(x, 4))
print("autovalori:", np.round(np.sort(np.linalg.eigvals(P).real), 4))
def tempo(da, a):                       # tempo medio per andare da 'da' ad 'a'
    idx = [i for i in range(N + 1) if i != a]
    A = np.eye(N) - P[np.ix_(idx, idx)]
    return np.linalg.solve(A, np.ones(N))[idx.index(da)]
print("0 -> 2:", tempo(0, 2), " 0 -> 4:", tempo(0, 4))
L = (np.eye(N + 1) + P) / 2
print("pigra, n=20:", np.round(np.eye(N + 1)[0] @ np.linalg.matrix_power(L, 20), 4))

Output (estratto):

pi P == pi: True   ritorni: [16.  4.  2.6667  4.  16.]
autovalori: [-1.  -0.5  0.  0.5  1.]
0 -> 2: 2.6667   0 -> 4: 21.3333
pigra, n=20: [0.0633 0.2516 0.375  0.2484 0.0617]