Algoritmo di Metropolis-Hastings: differenze tra le versioni

Versione delle 20:13, 19 feb 2019

Versione attuale delle 15:46, 29 apr 2025

(13 versioni intermedie di 10 utenti non mostrate)

Riga 1:

{{torna a|Catena di Markov Monte Carlo}}

L''''algoritmo di Metropolis-Hastings''' è un metodo [[MCMC]] usato per generare dei valori x<sub>1</submath>x_1, x₂x_2, ..\ldots, x<sub>nx_n</submath> che presentano una [[distribuzione di probabilità|distribuzione]] <math>p(x)</math> fissata a priori. Non necessita che la funzionedistribuzione <math>p(x)</math> sia nota, è sufficiente anzi che sia conosciuta una funzione <math>f(x)</math> proporzionale a <math>p(x).</math> 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, come è spesso il caso nell'[[inferenza bayesiana]].

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

Riga 7:

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 <math>p(x)</math> voluta. Necessita di una funzione <math>f(x) \propto p(x)</math> e di una ''proposal distribution'' <math>J(x^*|x_i)</math> simmetrica, che rispetti cioè la proprietà <math>J(x^*|x_i) = J(x_i|x^*)</math>. Le scelte più comuni per la distribuzione di proposta sono la normale <math>\mathcal N(x_i, \delta^2)</math>e l'uniforme <math>\text{unif}(x_i-\delta, x_i+\delta)</math>, dove <math>\delta</math> è un parametro da specificare prima della partenza dell'algoritmo.

Ciascuna iterazione dell'algoritmo di [[Nicholas Constantine Metropolis|Metropolis]] consiste nei seguenti passaggi:

Riga 14:

# si calcola il rapporto <math>w = \frac{f(x^*)}{f(x_i)} = \frac{p(x^*)}{p(x_i)}</math>;

# se <math>w \geq 1</math> si accetta il nuovo valore <math>x^* = x_{i+1}</math>;

# se invece <math>w < 1</math> il nuovo valore deve essere accettato con probabilità <math>w</math>. Si genera quindi un numero random <math>r</math> distribuito uniformemente nell'intervallo <math>[0, 1]</math>;

## se <math>r < w</math> si accetta il nuovo valore <math>x^* = x_{i+1}</math>;

## altrimenti il nuovo valore viene rigettato e si pone <math>x_{i+1} = x_i</math>.

Per generare una sequenza di ''<math>N''</math> elementi basta ripetere queste operazioni <math>N</math> volte a partire da un valore iniziale x<submath>0x_0,</submath>, scelto arbitrariamente.

Per avere una buona stima di <math>p(x)</math> è necessario generare sequenze abbastanza lunghe. La scelta del parametro δ<math>\delta</math> 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 [[autocorrelazione|autocorrelati]].

Di conseguenza, essendo δ<math>\delta</math> dipendente dalla forma e dalla scala di <math>p(x),</math> 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 <math>x</math> tali che <math>p(x)</math> 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 <math>x^*</math> e accettarlo o respingerlo in toto, è possibile considerare a parte ogni elemento di <math>x = (x₁x_1, ... x<sub>p\ldots,x_n)</submath>) e generare a parte un nuovo valore per ciascuno di questi elementi tramite una distribuzione simmetrica J<sub>j</submath>J_j(x^*_j_j|x<sub>jx_j)</submath>) per poi accettare o respingere questo valore singolarmente, al fine di definire x<submath>x_{i+1}.</submath>.

== 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 ''<math>w''</math>, che si ridefinisce come <math>w = \frac{J(x_i|x^*)}{J(x_i|x^*|x_i)}\frac{f(x^*)}{f(x_i)} = \frac{J(x_i|x^*)}{J(x_i|x^*|x_i)}\frac{p(x^*)}{p(x_i)}</math>. Il resto dell'algoritmo rimane invariato.

== Derivazione formale ==

Lo scopo dell'algoritmo di Metropolis-Hastings è generare una collezione di stati secondo una distribuzione desiderata <math>P(x)</math>. Per farlo, l'algoritmo utilizza un [[processo di Markov]], che asintoticamente raggiunge una [[Processo markoviano#Distribuzioni stazionarie|distribuzione stazionaria]] unica <math>\pi(x)</math> tale che <math>\pi(x) = P(x)</math>.

Un processo di Markov è definito univocamente dalle sue probabilità di transizione <math>P(x' \mid x)</math>, cioè la probabilità di passare da uno stato dato <math>x</math> a un altro stato dato <math>x'</math>. Esso ammette una distribuzione stazionaria unica <math>\pi(x)</math> quando sono soddisfatte le seguenti due condizioni:

# ''Esistenza della distribuzione stazionaria'': deve esistere una distribuzione stazionaria <math>\pi(x)</math>. Una condizione sufficiente (ma non necessaria) è il detailed balance (bilanciamento dettagliato), che richiede che ogni transizione <math>x \to x'</math> sia reversibile, cioè per ogni coppia di stati <math>x, x'</math>, la probabilità di essere nello stato <math>x</math> e passare a <math>x'</math> deve essere uguale alla probabilità di essere nello stato <math>x'</math> e passare a <math>x</math>, ovvero: <math>\pi(x)P(x'\mid x)=\pi(x')P(x\mid x')</math>. Una distribuzione <math>\pi(x)</math> che soddisfa il bilancio dettagliato è necessariamente stazionaria, infatti segue che:

::<math>\int \pi(x)P(x'\mid x) dx= \int \pi(x')P(x\mid x') dx = \pi(x')\int P(x\mid x') dx = \pi(x'),</math>

:cioè <math>\pi(x)</math> soddisfa la definizione di distribuzione stazionaria.

# ''Unicità della distribuzione stazionaria'': la distribuzione stazionaria <math>\pi(x)</math> deve essere unica. Questo è garantito dall'[[Processo markoviano#Catene di Markov ergodiche|ergodicità]] del processo di Markov, la quale richiede che ogni stato sia:

## ''aperiodico'' - il sistema non ritorna allo stesso stato a intervalli fissi;

## ''positivamente ricorrente'' - il numero atteso di passi per tornare nello stesso stato è finito.

L'algoritmo di Metropolis–Hastings consiste nel progettare un processo di Markov (costruendo le [[probabilità di transizione]]) che soddisfi le due condizioni sopra, in modo che la sua distribuzione stazionaria <math>\pi(x)</math> sia esattamente <math>P(x)</math>. La derivazione dell'algoritmo parte dalla condizione di bilanciamento dettagliato:

:<math>P(x)P(x'\mid x)=P(x')P(x\mid x')</math>

che può essere riscritta come:

:<math>\frac{P(x' \mid x)}{P(x \mid x')} = \frac{P(x')}{P(x)}.</math>

L'approccio consiste nel suddividere la transizione in due sotto-passi: la ''proposta'' e l'''accettazione/rifiuto''. La distribuzione di proposta <math>g(x' \mid x)</math> è la probabilità condizionata di proporre uno stato <math>x'</math> dato <math>x</math>, mentre la distribuzione di accettazione <math>A(x', x)</math> è la probabilità di accettare lo stato proposto <math>x'</math>. La probabilità di transizione può essere scritta come prodotto di queste:

:<math>P(x'\mid x) = g(x' \mid x) A(x', x).</math>

Inserendo questa relazione nell'equazione precedente si ottiene:

:<math>\frac{A(x', x)}{A(x, x')} = \frac{P(x')}{P(x)}\frac{g(x \mid x')}{g(x' \mid x)}.</math>

Il passo successivo nella derivazione è scegliere un rapporto di accettazione che soddisfi la condizione sopra. Una scelta comune è quella di Metropolis:

:<math>A(x',x) = \begin{cases}

\min\left(1, \frac{P(x')}{P(x)} \frac{g(x \mid x')}{g(x' \mid x)}\right) & \text{se } P(x)g(x' \mid x) \ne 0; \\

1 & \text{altrimenti}.

\end{cases}</math>

Per questo rapporto di accettazione <math>A</math> la condizione è soddisfatta. La scelta di <math>A</math> siffatta è giustificata da due punti:

* Muoversi verso zone di più alta densità di probabilità aumenta la probabilità di accettazione, infatti in tal caso <math>\frac{P(x')}{P(x)} > 1</math>.

* Il termine <math>\frac{g(x \mid x')}{g(x' \mid x)}</math> corregge l'eventuale asimmetria della probabilità di transizione, così da rispettare il bilancio dettagliato, infatti, consideriamo il caso in cui <math>P(x)g(x'\mid x) > P(x')g(x\mid x')</math>, cioè è più probabile andare da <math>x</math> in <math>x'</math> piuttosto che da <math>x'</math> in <math>x</math>, allora <math>A(x',x) < 1</math> mentre <math>A(x, x')=1</math> e quindi

::<math>P(x)P(x' \mid x) = P(x)g(x'\mid x)A(x', x) = P(x')g(x\mid x')\underbrace{A(x, x')}_{=1} = P(x')P(x\mid x').</math>

:Quindi <math>P(x)</math> è effettivamente la distribuzione limitre della catena di Markov.

L'algoritmo di Metropolis-Hastings può quindi essere scritto come segue:

# Inizializzazione

## scegliamo uno stato iniziale <math>x_0</math>;

## poniamo <math>t = 0</math>.

# Iterazione

## generiamo un candidato <math>x'</math> secondo <math>g(x' \mid x_t)</math>_;

## calcoliamo la probabilità di accettazione <math>A(x', x_t)</math>.

# Accettazione o rifiuto

## generiamo un [[numero casuale]] uniforme <math>u \in [0, 1]</math>;

## se <math>u \le A(x' , x_t)</math>, allora accettiamo il candidato e assegnamo: <math>x_{t+1} = x'</math>;

## se <math>u > A(x' , x_t)</math>, allora rifiutiamo il candidato e assegnamo: <math>x_{t+1} = x_t</math>.

# Passiamo alla prossima iteriazione, <math>t = t+1</math>.

A condizione che siano soddisfatte le ipotesi richieste, la distribuzione empirica degli stati salvati <math>x_0, \dots, x_T</math> tenderà a <math>P(x)</math>. Il numero di iterazioni <math>(T)</math> necessario per stimare efficacemente <math>P(x)</math> dipende da numerosi fattori, tra cui la relazione tra <math>P(x)</math> e la distribuzione di proposta, e la precisione desiderata nella stima. Per distribuzioni su spazi discreti degli stati, deve essere dell'ordine del tempo di autocorrelazione del processo di Markov.

È importante notare che, in un problema generale, non è chiaro quale distribuzione <math>g(x' \mid x)</math> si debba usare, né quante iterazioni siano necessarie per una buona stima: entrambi sono parametri liberi del metodo, che devono essere adattati al problema specifico.

== 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 ==

* {{Cita libro|nome=Hoff, Peter|cognome=D.|titolo=A first course in Bayesian statistical methods|url=https://www.worldcat.org/oclc/432708578|accesso=2018-12-28|data=2009|editore=Springer|OCLC=432708578|ISBN=9780387924076}}

* {{Cita pubblicazione|autore= [[Nicholas Metropolis]]|etal= sì|titolo=Equation of State Calculations by Fast Computing Machines |rivista=The Journal of Chemical Physics |url=https://aip.scitation.org/action/captchaChallenge?redirectUrl=https%3A%2F%2Faip.scitation.org%2Fdoi%2F10.1063%2F1.1699114&|doi= 10.1063/1.1699114|anno=1953 |accesso=2018-12-28}}

* {{Cita pubblicazione|nome=W. K.|cognome=Hastings|data=1970-04-01|titolo=Monte Carlo sampling methods using Markov chains and their applications|rivista=Biometrika|volume=57|numero=1|pp=97–10997-109|lingua=en|accesso=2018-12-28|doi=10.1093/biomet/57.1.97|url=https://academic.oup.com/biomet/article/57/1/97/284580}}

* Raftery, Adrian E., and Steven Lewis. "[https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=0daf54c4b59fd2c362de822de0ffdab84f49c6fd How Many Iterations in the Gibbs Sampler?]" ''In Bayesian Statistics 4''. 1992.

== Voci correlate ==