Monte Carlo a Catena di Markov

Una catena di Markov è un processo stocastico discreto \(X_n\) in cui il futuro dipende solo dal presente. Tutti gli \(X_n\)hanno gli stessi valori, l'insieme finito degli stati, con una distribuzione iniziale \(𝜋_0.\) Con le giuste proprietà, una catena di Markov converge a una distribuzione stazionaria unica, indipendentemente dal punto di partenza.

Inserisci la matrice di transizione \(𝑃_{𝑖,𝑗}=𝑃(𝑋_{𝑛+1}=𝑗∣𝑋_𝑛=𝑖):\)
[ [0.2, 0.5, 0.2, 0.1], [0.1, 0.6, 0.2, 0.1], [0.3, 0.3, 0.3, 0.1], [0.0, 0.0, 0.0, 1.0] ]
Stato iniziale:    Passi per blocco:
Esegui blocco

Nella stragrande maggioranza dei problemi reali, quando per generare numeri casuali che seguano specifiche distribuzioni di probabilità non è possibile usare l'inversione della funzione di ripartizione, quando lo spazio degli stati è continuo, \(p(x)\) è nota solo a meno di una costante (caso bayesiano) e la costante di normalizzazione è un integrale impossibile, oppure quando lo spazio discreto è enorme e non si possono enumerare tutti gli stati e le loro probabilità, si usano algoritmi di una classe di metodi Monte Carlo a Catena di Markov (MCMC).
Data una distribuzione di probabilità, è possibile costruire una catena di Markov la cui distribuzione degli elementi la approssima, ovvero la distribuzione di equilibrio della catena di Markov corrisponde alla distribuzione target.
Si può dire anche che MCMC costruisce una passeggiata casuale tra i valori andando più volentieri verso i valori più frequenti e tornando indietro più volentieri dai valori rari. Dopo un po' la percentuale di tempo passata su ogni valore è uguale alla sua probabiltà originale.

Uno di questi metodi è l'algoritmo Metropolis-Hastings.

Ad esempio supponiamo di voler campionare da una distribuzione nota per cui sia costoso calcolare numericamente l'integrale per normalizzare. Con MCMC esplorariamo lo spazio delle sequenze senza enumerarlo esplicitamente.
Partendo da un punto qualsiasi nello spazio dei valori, ad ogni passo:

1. si propone uno spostamento verso un nuovo punto (secondo una regola qualsiasi);
2. si accetta o rifiuta quel punto con una probabilità calcolata in modo da rispettare la bilancia dettagliata — la condizione che garantisce la convergenza alla distribuzione target

Dopo molti passi, i campioni raccolti si comportano come se provenissero dalla distribuzione target.

Più precisamente se \(p(x)\) è la distribuzione nota:

1. dato lo stato attuale \(x;\)
2. campiona un candidato \(x^*\) da una distribuzione di proposta \(q(x^*|x);\)
3. calcola il rapporto di accettazione: \[\alpha = \min\!\left(1,\ \frac{p(x^*)\, q(x\mid x^*)}{p(x)\, q(x^*\mid x)}\right);\]
4. accetta \(x^*\) con probabilità \(α\), altrimenti rimane in \(x.\)

Non serve conoscere la costante di normalizzazione di \(p(x),\) perché si cancella nel rapporto. Questo è il motivo per cui MCMC è così potente in statistica bayesiana, dove spesso quella costante (l'evidenza marginale) non è calcolabile.
I primi campioni dovranno essere scartati (Burn-in) poiché la catena non ha ancora raggiunto la stazionarietà. I campioni consecutivi sono correlati, non indipendenti.

La MCMC costruisce dunque una catena di Markov la cui distribuzione stazionaria è la distribuzione da campionare. La catena deve "esplorare" bene lo spazio — una proposta troppo piccola o troppo grande rallenta la convergenza. Si usano tracce, \(\hat R\) di Gelman-Rubin, e più catene parallele per verificare la convergenza.

Possiamo vedere come la distribuzione ottenuta con questo algoritmo approssimi la distribuzione da campionare.

Nel caso Metropolis-Hastings con proposta simmetrica e \(n\) stati la matrice di transizione \( P \), di dimensione \( n \times n \), dove l'elemento \( P_{ij} \) è la probabilità di passare dallo stato \( i \) allo stato \( j \) in un passo dell'algoritmo MCMC, si può descrivere come: \[ P_{ij} = Q_{ij} \cdot \alpha_{ij} \quad \text{per } j \neq i \] \[ P_{ii} = Q_{ii} + \sum_{k \neq i} Q_{ik}(1 - \alpha_{ik}) \] dove

• \( Q \) è la matrice delle proposte, con \( Q_{ij} \) probabilità di proporre \( j \) quando si è in \( i ,\) e può essere semplicemente \(Q_{ii}=0,\; Q_{ij}=\frac{1}{n} ,\)
• \( \alpha_{ij} = \min\left(1, \frac{p(j)}{p(i)}\right) \) è la probabilità di accettare.

Proprietà fondamentali sono:

• stazionarietà, per cui per la distribuzione target \( \pi \) normalizzata vale \[ \pi P = \pi; \]inoltre dopo \(t\) passi è \(π_t=π_0P^t\) e per \(t→ ∞\) si ha \(π_t→π;\)
• bilancio dettagliato (condizione sufficiente), cioè l'algoritmo è costruito per soddisfare \[π_iP_{ij}=π_jP_{ji},\;\; ∀i,j;\]
• la catena è irriducibile, cioè ogni stato è raggiungibile da ogni altro, e aperiodica, cioè si può restare nello stesso stato.

Varianti principali a Metropolis-Hastings sono:

• Gibbs Sampling, campiona ogni variabile condizionatamente alle altre, efficiente quando le condizionali sono note in forma chiusa;
• Hamiltonian Monte Carlo, usa la geometria del gradiente per fare proposte molto più efficienti ed è alla base di Stan e PyMC;
• No-U-Turn Sampler, versione adattiva dell'Hamiltonian Monte Carlo, oggi lo standard de facto.