Esercizio - Metropolis-Hastings su tre stati con proposta asimmetrica

Teoria: Metodi Monte Carlo - stima, generazione di variabili aleatorie e MCMCMonte Carlo stima una grandezza come media di campioni casuali: errore tipico σ/√N, numero di campioni N ≥ ln(2/δ)/(2ε²) per Hoeffding (probabilità), N ≥ 3 ln(2/δ)/(pε²) per Chernoff (errore relativo). Stima di π con p=π/4; integrali come media di g(U); eventi rari con campionamento per importanza (pesi φ/q); conteggio DNF con Karp-Luby in tempo polinomiale perché P(Y=1) ≥ 1/m. Generazione: trasformazione inversa X=F^{-1}(U) e accettazione-rifiuto (f ≤ M g, efficienza 1/M). MCMC: Metropolis-Hastings costruisce una catena con stazionaria π (nota a meno di costante) imponendo il bilancio dettagliato con probabilità di accettazione min(1, π(y)Q(y,x)/(π(x)Q(x,y))).Metodi Monte Carlo - stima, generazione di variabili aleatorie e MCMC →; 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. Si vuole campionare dalla distribuzione $\pi\propto w=(1,3,2)$ sugli stati $1,2,3$ (quindi $\pi=(\tfrac16,\tfrac12,\tfrac13)$). La proposta $Q$ è asimmetrica: da $1$ si propone sempre $2$, da $2$ si propone $1$ o $3$ con probabilità $\tfrac12$, da $3$ si propone sempre $2$.

1. Calcolare le probabilità di accettazione di Metropolis-Hastings e la matrice $P$ della catena.
2. Verificare il bilancio dettagliato e che $\pi P=\pi$.
3. Cosa si ottiene se per errore si usa la regola simmetrica $a=\min(1,w(y)/w(x))$ senza il fattore $Q(y,x)/Q(x,y)$?
4. Distanza di variazione totale da $\pi$ dopo $n$ passi partendo da $1$, per $n=1,2,3,5,10,15,20$; autovalori; quanti passi per portare la distanza sotto $0{,}01$?
5. Stimare $E_\pi[X]$ con la media lungo la catena.

Soluzione

1. La regola di Metropolis-Hastings è $a(x,y)=\min\Bigl(1,\dfrac{w(y)\,Q(y,x)}{w(x)\,Q(x,y)}\Bigr)$.

• $a(1,2)=\min\bigl(1,\tfrac{3\cdot Q(2,1)}{1\cdot Q(1,2)}\bigr)=\min\bigl(1,\tfrac{3\cdot1/2}{1\cdot1}\bigr)=\min(1,1{,}5)=1$.
• $a(2,1)=\min\bigl(1,\tfrac{1\cdot Q(1,2)}{3\cdot Q(2,1)}\bigr)=\min\bigl(1,\tfrac{1}{3\cdot1/2}\bigr)=\tfrac23$.
• $a(2,3)=\min\bigl(1,\tfrac{2\cdot Q(3,2)}{3\cdot Q(2,3)}\bigr)=\min\bigl(1,\tfrac{2\cdot1}{3\cdot1/2}\bigr)=\min(1,\tfrac43)=1$.
• $a(3,2)=\min\bigl(1,\tfrac{3\cdot Q(2,3)}{2\cdot Q(3,2)}\bigr)=\min\bigl(1,\tfrac{3\cdot1/2}{2\cdot1}\bigr)=\tfrac34$.

Le probabilità di transizione sono $P(x,y)=Q(x,y)a(x,y)$ per $y\ne x$, e il resto va sul cappio: $$P=\begin{pmatrix}0&1&0\\\tfrac13&\tfrac16&\tfrac12\\0&\tfrac34&\tfrac14\end{pmatrix}.$$ Per esempio $P(2,1)=\tfrac12\cdot\tfrac23=\tfrac13$, $P(2,3)=\tfrac12\cdot1=\tfrac12$, $P(2,2)=1-\tfrac13-\tfrac12=\tfrac16$, $P(3,2)=1\cdot\tfrac34$, $P(3,3)=\tfrac14$.

2. Bilancio dettagliato: $\pi_1P_{12}=\tfrac16\cdot1=\tfrac16$ e $\pi_2P_{21}=\tfrac12\cdot\tfrac13=\tfrac16$ ✓; $\pi_2P_{23}=\tfrac12\cdot\tfrac12=\tfrac14$ e $\pi_3P_{32}=\tfrac13\cdot\tfrac34=\tfrac14$ ✓; $\pi_1P_{13}=0=\pi_3P_{31}$ ✓. Quindi $\pi P=\pi$: per esempio la componente $3$ di $\pi P$ è $\tfrac12\cdot\tfrac12+\tfrac13\cdot\tfrac14=\tfrac14+\tfrac1{12}=\tfrac13$ ✓. La catena è irriducibile (si va $1\leftrightarrow2\leftrightarrow3$) e aperiodica (cappi in $2$ e $3$): converge a $\pi$ da ogni stato.

3. Con $a'(x,y)=\min(1,w(y)/w(x))$: $a'(1,2)=1$, $a'(2,1)=\tfrac13$, $a'(2,3)=\tfrac23$, $a'(3,2)=1$. Allora $P'(2,1)=\tfrac12\cdot\tfrac13=\tfrac16$, $P'(2,3)=\tfrac12\cdot\tfrac23=\tfrac13$, $P'(2,2)=\tfrac12$, $P'(3,2)=1$, $P'(3,3)=0$. Anche $P'$ è una catena nascita-morte e la sua stazionaria $\pi'$ si trova dal bilancio dettagliato: $\pi'_1\cdot1=\pi'_2\cdot\tfrac16$ e $\pi'_2\cdot\tfrac13=\pi'_3\cdot1$ danno $\pi'_1=\tfrac16\pi'_2$, $\pi'_3=\tfrac13\pi'_2$, somma $\tfrac32\pi'_2=1$, quindi $\pi'=(\tfrac19,\tfrac23,\tfrac29)\ne(\tfrac16,\tfrac12,\tfrac13)$. L'effetto della proposta non viene compensato e la media di $X$ risulta sbagliata ($\tfrac{1+12+6}{9}=2{,}111$ invece di $2{,}167$). Il fattore $Q(y,x)/Q(x,y)$ serve a compensare le asimmetrie della proposta.

4. Gli autovalori di $P$ sono $1$, $0{,}1187$ e $-0{,}7020$. Il secondo in modulo è $0{,}702$: l'errore decade come $0{,}702^n$ (con oscillazione di segno, perché è negativo). Partendo da $1$:

La distanza scende sotto $0{,}01$ tra $n=10$ ($0{,}0198$) e $n=15$ ($0{,}0034$): il primo $n$ con distanza $<0{,}01$ è $12$ ($0{,}0139$ a $n=11$, $0{,}0098$ a $n=12$). Per un uso pratico si scartano i primi $\approx12$–$20$ passi (burn-in).

5. Per $X$ in $\{1,2,3\}$: $E_\pi[X]=\tfrac16+2\cdot\tfrac12+3\cdot\tfrac13=\tfrac{13}6=2{,}1667$. Simulando la catena MH per $3\cdot10^5$ passi (codice sotto) si ottengono le frequenze $(0{,}1658;\ 0{,}4993;\ 0{,}3349)$ e la media $2{,}169$, in accordo.

Codice di verifica

import numpy as np
w = np.array([1., 3., 2.]); pi = w / w.sum()
Q = np.array([[0, 1, 0], [.5, 0, .5], [0, 1, 0]])

def matrice(regola):
    P = np.zeros((3, 3))
    for x in range(3):
        for y in range(3):
            if x != y and Q[x, y] > 0: P[x, y] = Q[x, y] * regola(x, y)
        P[x, x] = 1 - P[x].sum()
    return P

mh = lambda x, y: min(1, w[y] * Q[y, x] / (w[x] * Q[x, y]))
sbagliata = lambda x, y: min(1, w[y] / w[x])
P, Ps = matrice(mh), matrice(sbagliata)
print(P, "\npi P == pi:", np.allclose(pi @ P, pi), " autovalori:", np.sort(np.linalg.eigvals(P).real))
v = np.linalg.eig(Ps.T)[1][:, np.argmax(np.linalg.eig(Ps.T)[0].real)].real
print("stazionaria della regola sbagliata:", v / v.sum())
x0 = np.array([1., 0, 0])
for n in (1, 2, 3, 5, 10, 12, 15, 20):
    d = x0 @ np.linalg.matrix_power(P, n)
    print(n, np.round(d, 4), "TV =", round(0.5 * np.abs(d - pi).sum(), 5))

rng = np.random.default_rng(1); x, cont = 0, np.zeros(3)
for _ in range(300_000):
    y = rng.choice(3, p=Q[x])
    if rng.random() < mh(x, y): x = y
    cont[x] += 1
print(cont / 300_000, (cont / 300_000) @ [1, 2, 3], pi @ [1, 2, 3])