12 - Catene di Markov e processi stocastici

Oggi piove. Che tempo farà fra dieci giorni, e quanti giorni di pioggia ti aspetti in un anno? Se supponi che il tempo di domani dipenda solo da quello di oggi, hai già costruito una catena di Markov, e le due domande hanno una risposta che si calcola con una moltiplicazione di matrici e un sistema lineare. Questo articolo, il dodicesimo della serie Probabilità e Statistica, parte da un meteo a tre stati e arriva a tre usi diversi della stessa idea: il comportamento di lungo periodo, il gioco d'azzardo che finisce in rovina e il PageRank. È anche la base del campionatore Metropolis visto nell'articolo sulla stima bayesiana. Ogni risultato è calcolato in modo esatto e confermato con una simulazione a seme fisso.

Che cos'è una catena di Markov e come si scrive

Un processo stocastico è una successione di variabili casuali X₀, X₁, X₂, … che descrive un fenomeno nel tempo. Si dice catena di Markov (a tempo discreto e con un numero finito di stati) se la distribuzione del passo successivo, dato tutto il passato, dipende solo dallo stato attuale: P(X_t+1 = j | X_t = i, passato) = P(X_t+1 = j | X_t = i). Se questa probabilità non cambia con t, la catena è omogenea e si riassume in una matrice di transizione P, dove l'elemento della riga i e colonna j è P_ij, la probabilità di passare da i a j. Ogni riga è una distribuzione di probabilità: contiene numeri non negativi con somma 1.

Nel meteo i tre stati sono sole, nuvoloso, pioggia. Dopo un giorno di sole il sole torna con probabilità 0,7, il nuvoloso con 0,2, la pioggia con 0,1; le altre due righe sono scelte in modo plausibile ma arbitrario. Il listato mostra tre fatti. Il primo è l'equazione di Chapman-Kolmogorov: le probabilità a n passi sono le potenze della matrice, perché per andare da i a j in m + n passi si passa per uno stato intermedio, e la somma sugli stati intermedi è proprio il prodotto di matrici. Il secondo è che dopo una decina di passi le righe di Pⁿ sono quasi uguali fra loro: il punto di partenza è stato dimenticato. Il terzo è il valore a cui convergono, la distribuzione stazionaria π, che soddisfa π P = π e ha somma 1. Il programma la ricava risolvendo il sistema in frazioni esatte, così il risultato non dipende da arrotondamenti, e poi la confronta con la frequenza di visita di una lunga traiettoria simulata.

from fractions import Fraction as F
import numpy as np

stati = ["sole", "nuvoloso", "pioggia"]
P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])
assert np.allclose(P.sum(axis=1), 1)          # ogni riga e' una distribuzione

# Probabilita' a n passi: Chapman-Kolmogorov dice che P^(m+n) = P^m P^n
P2, P10 = np.linalg.matrix_power(P, 2), np.linalg.matrix_power(P, 10)
print("P^2 riga 'sole'     :", np.round(P2[0], 4))
print("P^10 riga 'sole'    :", np.round(P10[0], 4))
print("P^10 riga 'pioggia' :", np.round(P10[2], 4))

# Distribuzione stazionaria pi (pi P = pi, somma 1), in frazioni esatte:
# elimino l'ultima equazione (ridondante) e la sostituisco con la normalizzazione
M = [[F(str(P[j][i])) - (1 if i == j else 0) for j in range(3)] for i in range(3)]  # (P^T - I)
M[2] = [F(1), F(1), F(1)]
b = [F(0), F(0), F(1)]
for c in range(3):                              # eliminazione di Gauss-Jordan
    r = next(r for r in range(c, 3) if M[r][c] != 0)
    M[c], M[r], b[c], b[r] = M[r], M[c], b[r], b[c]
    piv = M[c][c]
    M[c] = [x / piv for x in M[c]]; b[c] /= piv
    for r in range(3):
        if r != c:
            f = M[r][c]
            M[r] = [x - f * y for x, y in zip(M[r], M[c])]; b[r] -= f * b[c]
print("pi esatta           :", [str(x) for x in b], "=", [round(float(x), 5) for x in b])
pi = np.array([float(x) for x in b])
print("controllo pi P = pi :", np.allclose(pi @ P, pi))

# Autovalori: 1 e due valori di modulo minore di 1
print("autovalori          :", np.round(np.sort(np.linalg.eigvals(P).real)[::-1], 4))

# Simulazione a seme fisso: frazione di giorni in ciascuno stato su una sola traiettoria
rng = np.random.default_rng(12)
N = 300_000
x = 0
visite = np.zeros(3)
cum = P.cumsum(axis=1)
u = rng.random(N)
for t in range(N):
    x = int(np.searchsorted(cum[x], u[t]))
    visite[x] += 1
print("frequenze simulate  :", np.round(visite / N, 5))
print("scarto massimo      :", f"{np.abs(visite / N - pi).max():.5f}")

# Output:
# P^2 riga 'sole'     : [0.57 0.25 0.18]
# P^10 riga 'sole'    : [0.4568 0.2825 0.2607]
# P^10 riga 'pioggia' : [0.4562 0.2827 0.2611]
# pi esatta           : ['21/46', '13/46', '6/23'] = [0.45652, 0.28261, 0.26087]
# controllo pi P = pi : True
# autovalori          : [1.     0.4732 0.1268]
# frequenze simulate  : [0.45557 0.28303 0.2614 ]
# scarto massimo      : 0.00095

Da una giornata di sole, dopo due giorni la probabilità di sole scende da 0,7 a 0,57, e dopo dieci giorni a 0,4568: la riga della pioggia dà 0,4562, quasi la stessa cosa. La soluzione esatta è π = (21/46, 13/46, 6/23) = (0,4565; 0,2826; 0,2609): nel lungo periodo il 45,7% dei giorni è soleggiato e il 26,1% piovoso, a prescindere da come è cominciata la serie. Su 300.000 giorni simulati lo scarto massimo dalle frequenze esatte è 0,00095, coerente con l'errore atteso per 300.000 passi correlati, che stimo intorno a 0,0015 (i giorni successivi sono correlati, quindi l'errore è maggiore di quello di un campione indipendente). Gli autovalori di P sono 1, 0,4732 e 0,1268: l'autovalore 1 corrisponde a π, e il secondo in ordine di grandezza regola la velocità con cui si dimentica il punto di partenza, come vedrai fra poco.

Quando la catena dimentica da dove è partita

Non tutte le catene hanno questo comportamento, e le condizioni sono due. Una catena è irriducibile se da ogni stato si può raggiungere ogni altro in un numero finito di passi: altrimenti l'insieme degli stati si spezza in parti che non comunicano. È aperiodica se il ritorno a uno stato non è vincolato a multipli di un intero maggiore di 1. Il teorema di convergenza per catene finite (Levin, Peres e Wilmer, Teorema 4.9) afferma che una catena irriducibile e aperiodica ha un'unica distribuzione stazionaria π e che Pⁿ(i, ·) tende a π per ogni stato di partenza i. Se la catena è solo irriducibile, π esiste ed è unica, ma la convergenza di Pⁿ può mancare.

L'ultimo blocco del listato seguente è il controesempio più piccolo: due stati A e B che si scambiano sempre (A passa a B con probabilità 1 e viceversa). La catena è irriducibile ma ha periodo 2, e le potenze di P non convergono: A^10 è l'identità, A^11 lo scambio. Ciò che converge è la media temporale: la media delle prime 1000 potenze vale 0,5 in ogni entrata, cioè la stazionaria (1/2, 1/2). Per questo nel meteo, e in generale quando usi una catena per simulare, è la frazione di tempo trascorso in ogni stato a convergere sempre, mentre la distribuzione al passo n converge solo in assenza di periodicità.

Quanto ci mette a dimenticare, e quanto ci mette a tornare

Il programma qui sotto misura due grandezze. La prima è il tempo medio di ritorno: partendo da uno stato, quanti passi servono in media per ripassarci. Per una catena irriducibile sullo spazio finito vale 1/π_i (è il lemma di Kac, nel testo di Levin, Peres e Wilmer): se il 45,7% dei giorni è soleggiato, tra un giorno di sole e il successivo passano in media 1/0,4565 = 2,19 giorni. La seconda è la distanza in variazione totale dalla stazionaria, TV(n), cioè la metà della somma dei valori assoluti delle differenze fra Pⁿ(i, ·) e π, presa sul peggiore stato di partenza: misura quanto una distribuzione al passo n è ancora distinguibile da π.

import numpy as np

P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])
pi = np.array([21/46, 13/46, 6/23])             # stazionaria, calcolata nel listato precedente

# 1) Tempo medio di ritorno: per una catena irriducibile vale 1/pi_i (formula di Kac)
rng = np.random.default_rng(5)
cum = P.cumsum(axis=1)
def tempo_ritorno(stato):
    x, t = stato, 0
    while True:
        x = int(np.searchsorted(cum[x], rng.random())); t += 1
        if x == stato:
            return t
for i, nome in enumerate(["sole", "nuvoloso", "pioggia"]):
    r = np.array([tempo_ritorno(i) for _ in range(40_000)])
    print(f"ritorno a {nome:9s}: simulato {r.mean():.3f} +- {r.std() / np.sqrt(len(r)):.3f}   1/pi = {1/pi[i]:.3f}")

# 2) Quanto in fretta si dimentica il punto di partenza: distanza in variazione totale
#    TV(n) = max sulle righe di (1/2) * somma |P^n[i, j] - pi[j]|
lam2 = np.sort(np.abs(np.linalg.eigvals(P)))[-2]
print(f"secondo autovalore in modulo: {lam2:.4f}")
Pn = np.eye(3)
for n in range(1, 9):
    Pn = Pn @ P
    tv = 0.5 * np.abs(Pn - pi).sum(axis=1).max()
    print(f"n={n}  TV={tv:.6f}  TV/lam2^n={tv / lam2**n:.4f}")

# 3) Bilancio dettagliato: pi_i P_ij == pi_j P_ji ?
flusso = pi[:, None] * P
print("catena del meteo reversibile:", np.allclose(flusso, flusso.T))
Q = np.array([[0.5, 0.5, 0.0, 0.0],             # passeggiata su 4 nodi in fila, con code ferme
              [0.25, 0.5, 0.25, 0.0],
              [0.0, 0.25, 0.5, 0.25],
              [0.0, 0.0, 0.5, 0.5]])
w, v = np.linalg.eig(Q.T)
q = np.real(v[:, np.argmax(np.real(w))]); q /= q.sum()
print("stazionaria di Q:", np.round(q, 4), " reversibile:", np.allclose(q[:, None] * Q, (q[:, None] * Q).T))

# 4) Periodicita': due stati che si alternano
A = np.array([[0.0, 1.0], [1.0, 0.0]])
print("A^10 =", np.linalg.matrix_power(A, 10).tolist(), " A^11 =", np.linalg.matrix_power(A, 11).tolist())
print("media delle prime 1000 potenze:", np.round(sum(np.linalg.matrix_power(A, k) for k in range(1, 1001)) / 1000, 3).tolist())

# Output:
# ritorno a sole     : simulato 2.194 +- 0.013   1/pi = 2.190
# ritorno a nuvoloso : simulato 3.537 +- 0.018   1/pi = 3.538
# ritorno a pioggia  : simulato 3.869 +- 0.024   1/pi = 3.833
# secondo autovalore in modulo: 0.4732
# n=1  TV=0.256522  TV/lam2^n=0.5421
# n=2  TV=0.126522  TV/lam2^n=0.5650
# n=3  TV=0.060522  TV/lam2^n=0.5712
# n=4  TV=0.028722  TV/lam2^n=0.5728
# n=5  TV=0.013602  TV/lam2^n=0.5733
# n=6  TV=0.006438  TV/lam2^n=0.5734
# n=7  TV=0.003047  TV/lam2^n=0.5734
# n=8  TV=0.001442  TV/lam2^n=0.5734
# catena del meteo reversibile: False
# stazionaria di Q: [0.1667 0.3333 0.3333 0.1667]  reversibile: True
# A^10 = [[1.0, 0.0], [0.0, 1.0]]  A^11 = [[0.0, 1.0], [1.0, 0.0]]
# media delle prime 1000 potenze: [[0.5, 0.5], [0.5, 0.5]]

I tempi di ritorno simulati (2,194 ± 0,013; 3,537 ± 0,018; 3,869 ± 0,024, con 40.000 ritorni ciascuno) coincidono con 1/π entro due errori standard: 2,190, 3,538, 3,833. La distanza TV scende da 0,2565 al primo passo a 0,0014 all'ottavo, e il rapporto TV(n)/|λ₂|ⁿ, con λ₂ = 0,4732 secondo autovalore in modulo, si stabilizza a 0,5734: il decadimento è geometrico e la sua ragione è il secondo autovalore. Quanti passi servono perché TV scenda sotto 0,01 (una definizione usuale di tempo di mixing)? Bastano 6, e il numero cresce in modo logaritmico se chiedi una precisione maggiore. Per catene con stati numerosi, il secondo autovalore vicino a 1 segnala una convergenza lenta: è il motivo per cui i campionatori MCMC con un passo troppo piccolo restano a lungo vicino al punto di partenza.

Bilancio dettagliato: la catena che si può invertire nel tempo

Una catena con distribuzione stazionaria π è reversibile se vale il bilancio dettagliato: π_i P_ij = π_j P_ji per ogni coppia di stati, cioè il flusso medio da i a j uguaglia quello da j a i. La condizione implica π P = π ed è la proprietà che si costruisce di proposito quando si progetta un campionatore: nell'algoritmo di Metropolis la probabilità di accettazione è scelta proprio per far valere il bilancio dettagliato rispetto alla distribuzione bersaglio (Levin, Peres e Wilmer, cap. 3). Il terzo blocco del listato precedente la verifica su due catene. Il meteo non è reversibile: π_sole P_{sole→nuvoloso} = 0,0913 contro π_nuvoloso P_{nuvoloso→sole} = 0,0848. La passeggiata su quattro nodi in fila, invece, lo è, con stazionaria (1/6, 1/3, 1/3, 1/6): ogni catena nascita-morte, nella quale si può muovere solo verso i vicini, lo è.

Questo spiega un'asimmetria che sorprende i principianti: una catena non reversibile può convergere benissimo, come il meteo. La reversibilità non è necessaria per avere una stazionaria, è solo una condizione sufficiente e comoda da controllare, perché verificarla richiede di guardare le coppie di stati una alla volta, invece di risolvere un sistema di equazioni.

Catene con stati assorbenti: la rovina del giocatore

Uno stato è assorbente se una volta raggiunto non si lascia più. In una catena con stati assorbenti conviene riordinare la matrice in blocchi: Q contiene le transizioni fra gli stati transitori e R quelle verso gli stati assorbenti. La matrice fondamentale N = (I − Q)⁻¹ ha come elemento (i, j) il numero atteso di visite allo stato transitorio j partendo da i. La somma di una riga di N è il numero atteso di passi prima dell'assorbimento, e N R dà la probabilità di finire in ciascuno stato assorbente. Il caso classico è la rovina del giocatore: parti con 3 euro, a ogni giro vinci un euro con probabilità p o ne perdi uno con probabilità 1 − p, e smetti a 0 o a 10.

import numpy as np

# Rovina del giocatore: capitale 0..N, a ogni giro +1 con prob. p, -1 con prob. 1-p; 0 e N assorbono
def rovina(N, p):
    n_trans = N - 1                              # stati transitori: 1..N-1
    Q = np.zeros((n_trans, n_trans)); R = np.zeros((n_trans, 2))   # R: [0, N]
    for i in range(1, N):
        if i + 1 < N: Q[i-1, i] = p
        else:         R[i-1, 1] = p
        if i - 1 > 0: Q[i-1, i-2] = 1 - p
        else:         R[i-1, 0] = 1 - p
    Nf = np.linalg.inv(np.eye(n_trans) - Q)      # matrice fondamentale: visite attese a ogni stato
    return Nf @ R, Nf.sum(axis=1)                # probabilita' di assorbimento, durata attesa

N = 10
for p in (0.5, 0.45):
    ass, durata = rovina(N, p)
    print(f"p={p}: da 3 euro P(arrivo a {N}) = {ass[2, 1]:.5f}, durata attesa = {durata[2]:.4f} giri")
print("p=0.5, formule chiuse: i/N =", 3/N, " e i(N-i) =", 3 * (N - 3))

rng = np.random.default_rng(31)
def partita(i, N, p):
    t = 0
    while 0 < i < N:
        i += 1 if rng.random() < p else -1; t += 1
    return i == N, t
R_ = np.array([partita(3, N, 0.45) for _ in range(100_000)])
v = R_[:, 0].mean(); d = R_[:, 1].astype(float)
print(f"p=0.45 simulato: P(arrivo) = {v:.5f} +- {np.sqrt(v*(1-v)/len(R_)):.5f}, durata = {d.mean():.3f} +- {d.std()/np.sqrt(len(d)):.3f}")

# Output:
# p=0.5: da 3 euro P(arrivo a 10) = 0.30000, durata attesa = 21.0000 giri
# p=0.45: da 3 euro P(arrivo a 10) = 0.12825, durata attesa = 17.1748 giri
# p=0.5, formule chiuse: i/N = 0.3  e i(N-i) = 21
# p=0.45 simulato: P(arrivo) = 0.12860 +- 0.00106, durata = 17.261 +- 0.055

Con una moneta onesta (p = 0,5) il calcolo matriciale restituisce 0,30000 di probabilità di arrivare a 10 e una durata attesa di 21 giri, esattamente le formule chiuse i/N e i(N − i) (Levin, Peres e Wilmer, §2.1). Con p = 0,45, un vantaggio del banco di soli cinque punti percentuali sotto il 50%, la probabilità di arrivare a 10 crolla a 0,12825 e la durata attesa scende a 17,17 giri. La simulazione con 100.000 partite lo conferma: 0,12860 ± 0,00106 e 17,261 ± 0,055. Il messaggio pratico è che una piccola asimmetria per giro diventa una grande asimmetria sul risultato finale, ed è la ragione per cui insistere a giocare contro un banco con un lieve vantaggio finisce quasi sempre in rovina.

PageRank: una catena per ordinare le pagine

Il PageRank (Brin e Page, 1998) assegna a ogni pagina web l'importanza che avrebbe un navigatore casuale: da una pagina segue uno dei suoi link con probabilità d (0,85 nel lavoro originale), oppure salta a una pagina qualsiasi con probabilità 1 − d. L'importanza è la distribuzione stazionaria di questa catena. Il salto casuale serve a una ragione matematica precisa: garantisce che la catena sia irriducibile e aperiodica anche se il grafo del web ha pagine isolate o gruppi di pagine che non si linkano, e quindi che π esista, sia unica e si possa calcolare iterando il prodotto x ← x M (il metodo delle potenze) da qualunque punto di partenza. Il listato lo fa su cinque pagine e verifica il risultato con l'autovettore sinistro.

import numpy as np

# PageRank su 5 pagine: link[i] = pagine verso cui punta i
link = {0: [1, 2], 1: [2], 2: [0], 3: [2, 4], 4: [3]}
n, d_ = 5, 0.85                                   # d_: probabilita' di seguire un link
G = np.zeros((n, n))
for i, out in link.items():
    for j in out: G[i, j] = 1 / len(out)
M = d_ * G + (1 - d_) / n                         # con prob. 1-d_ si salta a caso: la catena e' irriducibile e aperiodica
x = np.full(n, 1 / n)
for it in range(1, 200):
    nuovo = x @ M
    if np.abs(nuovo - x).sum() < 1e-12: break
    x = nuovo
w, v = np.linalg.eig(M.T)
e = np.real(v[:, np.argmax(np.real(w))]); e /= e.sum()
print(f"iterazione di potenza: {it} passi", np.round(x, 5))
print("autovettore sinistro  :", np.round(e, 5), " uguali:", np.allclose(x, e))

# Output:
# iterazione di potenza: 55 passi [0.32731 0.16911 0.34977 0.08689 0.06693]
# autovettore sinistro  : [0.32731 0.16911 0.34977 0.08689 0.06693]  uguali: True

L'iterazione converge in 55 passi, con tolleranza 10^-12 sulla somma delle differenze, e coincide con l'autovettore: la pagina 2, linkata da 0, 1 e 3, ha l'importanza maggiore (0,3498), seguita da 0 (0,3273). Le pagine 3 e 4, che si linkano soltanto fra loro e non ricevono nessun link dalle altre tre, restano in fondo (0,0869 e 0,0669). Il numero di passi dipende da d: più d è vicino a 1 più la convergenza è lenta, perché cala il peso del salto casuale che rimescola. Questo è un modello del PageRank, non il sistema di ranking di un motore di ricerca, che usa molti altri segnali e non è descritto da questo articolo.

Esercizio svolto e prova da fare da solo

Una macchina ha due stati, funzionante (F) e guasta (G). Ogni giorno una macchina funzionante si guasta con probabilità 0,1 e una guasta viene riparata con probabilità 0,4. La stazionaria soddisfa π_F · 0,1 = π_G · 0,4 (bilancio dettagliato, valido per ogni catena a due stati): π_G = π_F / 4, quindi π_F = 0,8 e π_G = 0,2. La macchina è guasta un giorno su cinque. Il tempo medio di ritorno al guasto è 1/0,2 = 5 giorni, mentre ogni guasto dura in media 1/0,4 = 2,5 giorni.

Che cosa sai fare adesso

Sai scrivere la matrice di transizione di un fenomeno, calcolare le probabilità a n passi con le potenze, ricavare la distribuzione stazionaria in forma esatta e riconoscere quando irriducibilità e aperiodicità ne garantiscono la convergenza. Sai usare il tempo di ritorno 1/π e il secondo autovalore per stimare la velocità di mixing, controllare il bilancio dettagliato e risolvere una catena con stati assorbenti con la matrice fondamentale. Questi strumenti compaiono nei campionatori MCMC dell'articolo precedente, e nel processo decisionale di Markov del reinforcement learning. Per il lato machine learning vedi anche Catene di Markov: la proprietà che dimentica il passato. Il prossimo articolo della serie passa ai test non parametrici.

Nella stessa serie trovi anche Analisi dati esplorativa (EDA) con Python.