Algoritmo di Metropolis-Hastings: differenze tra le versioni

Versione delle 00:21, 31 ago 2024

Versione attuale delle 15:46, 29 apr 2025

(7 versioni intermedie di 6 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 <math>x_1, x_2,\ldots, x_n</math> che presentano una [[distribuzione di probabilità|distribuzione]] <math>p(x)</math> fissata a priori. Non necessita che la distribuzione <math>p(x)</math> sia nota, è sufficiente 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 altrimentiin 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 32:

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^*|x_i)}\frac{f(x^*)}{f(x_i)} = \frac{J(x_i|x^*)}{J(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 ==

Riga 39 ⟶ 103:

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