Algoritmo di Metropolis-Hastings

L'algoritmo di Metropolis-Hastings è un metodo MCMC usato per generare dei valori x₁, x₂, .., x_n che presentano una distribuzione p(x) fissata a priori. Non necessita che la funzione p(x) sia nota, è sufficiente anzi che sia conosciuta una funzione f(x) proporzionale a p(x). Questo requisito così debole permette di usare l'algoritmo di Metropolis-Hastings per campionare da distribuzioni di cui l'integrale sia troppo difficile, o impossibile, da calcolare in forma analitica.

Il metodo è stato descritto da Hastings nel 1970, come generalizzazione dell'algoritmo di Metropolis del 1953.

Algoritmo di Metropolis

Per comprendere l'algoritmo generale è utile imparare prima quello originale, detto di Metropolis.

Il metodo si basa sulla generazione di valori "proposti" che vengono accettati o rigettati in modo da convergere alla distribuzione $p(x)$ voluta. Necessita di una funzione $f(x)\propto p(x)$ e di una proposal distribution $J(x^{*}|x_{i})$ simmetrica, che rispetti cioè la proprietà $J(x^{*}|x_{i})=J(x_{i}|x^{*})$. Le scelte più comuni per la distribuzione di proposta sono la normale ${\mathcal {N}}(x_{i},\delta ^{2})$e l'uniforme $unif(x_{i}-\delta ,x_{i}+\delta )$, dove delta è un parametro da specificare prima della partenza dell'algoritmo.

Ciascuna iterazione dell'algoritmo di Metropolis consiste nei seguenti passaggi:

1. si estrae un nuovo valore $x^{*}$ dalla distribuzione di proposta $J(x^{*}|x_{i})$;
2. si calcola il rapporto $w={\frac {f(x^{*})}{f(x_{i})}}={\frac {p(x^{*})}{p(x_{i})}}$;
3. se $w\geq 1$ si accetta il nuovo valore $x^{*}=x_{i+1}$;
4. se invece $w<1$ il nuovo valore deve essere accettato con probabilità $w$. Si genera quindi un numero random $r$ distribuito uniformemente nell'intervallo [0, 1];
   1. se $r<w$ si accetta il nuovo valore $x^{*}=x_{i+1}$;
   2. altrimenti il nuovo valore viene rigettato e si pone $x_{i+1}=x_{i}$.

Per generare una sequenza di N elementi basta ripetere queste operazioni N volte a partire da un valore iniziale x₀, scelto arbitrariamente.

Per avere una buona stima di p(x) è necessario generare sequenze abbastanza lunghe. La scelta del parametro δ può essere cruciale, se è troppo grande solo una piccola parte dei valori di prova proposti verrà accettato. Se invece è troppo piccolo la catena si muoverà molto lentamente e i valori risulteranno estremamente autocorrelati.

Di conseguenza, essendo δ dipendente dalla forma e dalla scala di p(x), deve essere di volta in volta calibrato correttamente; per la sua stima si può procedere per approssimazione successiva in modo che, fissato un delta, il numero di valori accettati sia un terzo del totale. Anche la scelta del valore iniziale è molto importante, in genere conviene partire da valori di x tali che p(x) assuma valori massimi in modo da avere una buona statistica nelle zone più probabili.

Caso multivariato

L'algoritmo descritto sopra funziona esattamente sia nel caso uni- che multivariato, ma esiste un secondo approccio al caso multivariato, particolarmente interessante quando si va a studiare la generalizzazione di Metropolis-Hastings. Anziché generare ad ogni iterazione un nuovo vettore x* e accettarlo o respingerlo in toto, è possibile considerare a parte ogni elemento di x = (x₁, ... x_p) e generare a parte un nuovo valore per ciascuno di questi elementi tramite una distribuzione simmetrica J_j(x*_j|x_j) per poi accettare o respingere questo valore singolarmente, al fine di definire x_i+1.

Algoritmo di Metropolis-Hastings

L'algoritmo di Metropolis richiede, per garantirne la convergenza limite, che la distribuzione di proposta sia simmetrica. Questa condizione limita di fatto il processo che genera i valori proposti al dominio dei random walk. Hastings (1970) propose una generalizzazione dell'algoritmo di Metropolis che permette la scelta di qualsiasi tipo di proposta.

L'algoritmo di Metropolis-Hastings procede nello stesso modo del suo predecessore, ma non richiede la simmetria della proposal distribution. Questo rilassamento delle ipotesi richiede una modifica nella definizione del rapporto w, che si ridefinisce come $w={\frac {J(x_{i}|x^{*})}{J(x_{i}|x^{*})}}{\frac {f(x^{*})}{f(x_{i})}}={\frac {J(x_{i}|x^{*})}{J(x_{i}|x^{*})}}{\frac {p(x^{*})}{p(x_{i})}}$. Il resto dell'algoritmo rimane invariato.

Tempi caratteristici

Affinché l'algoritmo perda memoria del dato iniziale e converga verso la distribuzione che si vuole campionare, è necessario eseguire un numero iniziale di iterazioni: tale numero si definisce tempo di termalizzazione. Similmente, nel calcolo degli errori è necessario considerare un tempo di correlazione, che consideri l'autocorrelazione tra due campionamenti successivi.

Bibliografia

• Hoff, Peter D., A first course in Bayesian statistical methods, Springer, 2009, ISBN 9780387924076, OCLC 432708578. URL consultato il 28 dicembre 2018.
• Nicholas Metropolis et al., Equation of State Calculations by Fast Computing Machines, in The Journal of Chemical Physics, 1953, DOI:10.1063/1.1699114. URL consultato il 28 dicembre 2018.
• (EN) W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, in Biometrika, vol. 57, n. 1, 1º aprile 1970, pp. 97–109, DOI:10.1093/biomet/57.1.97. URL consultato il 28 dicembre 2018.

Voci correlate

• Processo markoviano
• Nicholas Constantine Metropolis