Come funziona Hamiltonian Monte Carlo?
Ho realizzato il grafico seguente per spiegare come attualmente comprendo l'algoritmo HMC. Vorrei una verifica da un esperto in materia se questa comprensione è o non è corretta. Il testo nella diapositiva sottostante viene copiato di seguito per facilità di accesso:
Hamiltoniano Monte Carlo: un satellite orbita attorno a un pianeta. Più il satellite è vicino al pianeta, maggiori saranno gli effetti della gravità. Ciò significa, (A) maggiore energia potenziale e (B) maggiore energia cinetica necessaria per sostenere l'orbita. Quella stessa energia cinetica a una distanza maggiore dal pianeta, espellerebbe il satellite dall'orbita. Il satellite ha il compito di raccogliere foto di una specifica regione geografica. Più il satellite orbita intorno al pianeta, più velocemente si muove in orbita, più volte passa sopra la regione, più fotografie raccoglie. Al contrario, più un satellite è lontano dal pianeta, più lentamente si muove in orbita, meno volte passa sopra la regione, meno fotografie raccoglie. Nel contesto del campionamento, la distanza dal pianeta rappresenta la distanza dall'aspettativa di distribuzione. Un'area di bassa probabilità è lontana dalle aspettative; quando "orbita attorno a questa probabilità", un'energia cinetica inferiore significa meno campioni raccolti in un intervallo di tempo fisso, mentre quando orbita una probabilità maggiore significa più campioni raccolti dato lo stesso intervallo di tempo fisso. In una data orbita, l'energia totale, cinetica e potenziale, è costante; tuttavia, il rapporto tra i due non è semplice. Le equazioni hamiltoniane mettono in relazione i cambiamenti l'una nell'altra. Vale a dire, il gradiente di posizione rispetto al tempo è uguale alla quantità di moto. E il gradiente di quantità di moto rispetto al tempo è uguale al gradiente di energia potenziale rispetto alla posizione. Per calcolare la distanza percorsa da un satellite lungo il suo percorso orbitale, è necessario utilizzare l'integrazione del cavalletto, aggiornando in modo iterativo la quantità di moto e i vettori di posizione. Nel contesto del campionamento, la probabilità è analoga alla distanza dal pianeta e il gradiente di energia potenziale rispetto alla posizione è il gradiente della funzione di densità di probabilità rispetto al suo parametro di input, x. Questa informazione consente di esplorare il percorso orbitale attorno a vari input, X, corrispondenti alla stessa probabilità, y.
Tuttavia, non siamo semplicemente interessati ad esplorare una probabilità, dobbiamo esplorare più percorsi orbitali. Per ottenere ciò, la quantità di moto deve essere aumentata in modo casuale, portando il satellite più vicino o più lontano dal pianeta. Questi "calci di slancio" casuali consentono diverse probabilità di orbitare. Fortunatamente, le equazioni hamiltoniane assicurano che, indipendentemente dalla probabilità, il numero di campioni raccolti sia proporzionale alla probabilità, quindi i campioni raccolti seguono la forma della distribuzione target.
La mia domanda è: è questo un modo accurato di pensare a come funziona l'Hamiltoniano Monte Carlo?
Modificare:
L'ho implementato in un codice basato sulla mia comprensione dell'algoritmo. Funziona per un gaussiano con mu = 0, sigma = 1. Ma se cambio sigma si rompe. Ogni approfondimento è apprezzato.
import numpy as np
import random
import scipy.stats as st
import matplotlib.pyplot as plt
from autograd import grad
def normal(x,mu,sigma):
numerator = np.exp((-(x-mu)**2)/(2*sigma**2))
denominator = sigma * np.sqrt(2*np.pi)
return numerator/denominator
def neg_log_prob(x,mu,sigma):
num = np.exp(-1*((x-mu)**2)/2*sigma**2)
den = sigma*np.sqrt(np.pi*2)
return -1*np.log(num/den)
def HMC(mu=0.0,sigma=1.0,path_len=1,step_size=0.25,initial_position=0.0,epochs=1_000):
# setup
steps = int(path_len/step_size) -1 # path_len and step_size are tricky parameters to tune...
samples = [initial_position]
momentum_dist = st.norm(0, 1)
# generate samples
for e in range(epochs):
q0 = np.copy(samples[-1])
q1 = np.copy(q0)
p0 = momentum_dist.rvs()
p1 = np.copy(p0)
dVdQ = -1*(q0-mu)/(sigma**2) # gradient of PDF wrt position (q0) aka momentum wrt position
# leapfrog integration begin
for s in range(steps):
p1 += step_size*dVdQ/2 # as potential energy increases, kinetic energy decreases
q1 += step_size*p1 # position increases as function of momentum
p1 += step_size*dVdQ/2 # second half "leapfrog" update to momentum
# leapfrog integration end
p1 = -1*p1 #flip momentum for reversibility
#metropolis acceptance
q0_nlp = neg_log_prob(x=q0,mu=mu,sigma=sigma)
q1_nlp = neg_log_prob(x=q1,mu=mu,sigma=sigma)
p0_nlp = neg_log_prob(x=p0,mu=0,sigma=1)
p1_nlp = neg_log_prob(x=p1,mu=0,sigma=1)
# Account for negatives AND log(probabiltiies)...
target = q0_nlp - q1_nlp # P(q1)/P(q0)
adjustment = p1_nlp - p0_nlp # P(p1)/P(p0)
acceptance = target + adjustment
event = np.log(random.uniform(0,1))
if event <= acceptance:
samples.append(q1)
else:
samples.append(q0)
return samples
Ora funziona qui:
mu, sigma = 0,1
trial = HMC(mu=mu,sigma=sigma,path_len=2,step_size=0.25)
# What the dist should looks like
lines = np.linspace(-6,6,10_000)
normal_curve = [normal(x=l,mu=mu,sigma=sigma) for l in lines]
# Visualize
plt.plot(lines,normal_curve)
plt.hist(trial,density=True,bins=20)
plt.show()
Ma si rompe quando cambio sigma a 2.
# Generate samples
mu, sigma = 0,2
trial = HMC(mu=mu,sigma=sigma,path_len=2,step_size=0.25)
# What the dist should looks like
lines = np.linspace(-6,6,10_000)
normal_curve = [normal(x=l,mu=mu,sigma=sigma) for l in lines]
# Visualize
plt.plot(lines,normal_curve)
plt.hist(trial,density=True,bins=20)
plt.show()
Qualche idea? Mi sento come se fossi vicino a "ottenerlo".
Risposte
Prima di rispondere alla domanda su un modo intuitivo di pensare all'Hamiltoniano di Monte Carlo, è probabilmente meglio avere una solida conoscenza del normale MCMC. Mettiamo da parte la metafora del satellite per ora.
MCMC è utile quando desideri un campione non dichiarato da una distribuzione in cui hai solo qualcosa disponibile che è proporzionale al PDF, ma non il PDF stesso. Ciò si verifica (ad esempio) nelle simulazioni fisiche: il PDF è dato dalla distribuzione Boltzmann, p ~ exp (-E / kT), ma la cosa che puoi calcolare per qualsiasi configurazione del sistema è E, non p. La costante di proporzionalità non è nota, perché l'integrale di exp (-E / kT) sull'intero spazio di configurazione possibile è solitamente troppo difficile da calcolare. MCMC risolve questo problema facendo una passeggiata aleatoria in un modo specifico, dove la probabilità di prendere ("accettare") ogni passo è correlata al rapporto dei valori p (la costante di proporzionalità annulla). Nel tempo, la distribuzione dei campioni accettati dalla passeggiata aleatoria converge al PDF che vogliamo, senza mai dover calcolare esplicitamente p.
Nota che in quanto sopra, qualsiasi metodo per eseguire passi casuali è ugualmente valido, purché il camminatore casuale possa esplorare l'intero spazio. Il criterio di accettazione garantisce che i campioni selezionati convergono al PDF reale. In pratica, viene utilizzata una distribuzione gaussiana attorno al campione corrente (e il sigma può essere variato in modo che la frazione di passaggi accettati rimanga relativamente alta). In linea di principio non ci sarebbe nulla di sbagliato nel prendere passi da qualsiasi altra distribuzione continua ("distribuzione saltante") attorno al campione corrente, sebbene la convergenza possa essere molto più lenta.
Ora, Hamiltoniano Monte Carlo estende la metafora della fisica cercando specificamente di compiere passi in una direzione che è più probabile che venga accettata di un passo gaussiano. I passaggi sono ciò che farebbe un integratore a balestra, se cercasse di risolvere il movimento di un sistema in cui l'energia potenziale era E. Queste equazioni del moto includono anche un termine di energia cinetica, con una "massa" (non letteralmente fisica) e "quantità di moto". I passaggi che l'integratore leapfrog esegue in "tempo" vengono quindi passati come proposte all'algoritmo MCMC.
Perché funziona? La MC gaussiana esegue passi della stessa distanza in ogni direzione con uguale probabilità; l'unica cosa che lo spinge verso aree più densamente popolate del PDF è che è più probabile che i passi nella direzione sbagliata vengano rifiutati. L'Hamiltoniano MC propone passi sia nella direzione del gradiente E, sia nella direzione del moto accumulato nei passi recenti (direzione e grandezza della "quantità di moto"). Ciò consente un'esplorazione più rapida dello spazio e anche una maggiore probabilità di raggiungere più rapidamente regioni più densamente popolate.
Ora, la metafora del satellite: penso che questo non sia un modo molto utile per pensarci. I satelliti si muovono in un'orbita esatta; quello che hai qui è abbastanza casuale, più simile a una particella di gas in un contenitore con altre particelle. Ogni collisione casuale ti dà un "passo"; nel tempo la particella sarà ovunque nel contenitore con uguale probabilità (poiché il PDF qui è uguale ovunque, tranne le pareti che rappresentano un'energia molto alta / effettivamente zero PDF). L'MCMC gaussiano è come una particella di massa zero che fa una passeggiata casuale (o particella di massa diversa da zero in un mezzo relativamente viscoso): ci arriverà attraverso il moto browniano, ma non necessariamente veloce. L'Hamiltoniano MC è una particella con una massa diversa da zero: può raccogliere abbastanza quantità di moto per continuare ad andare nella stessa direzione nonostante le collisioni, e quindi a volte può sparare da un'estremità all'altra del contenitore (a seconda della sua massa rispetto alla frequenza / entità delle collisioni). Ovviamente rimbalzerebbe ancora contro i muri, ma in generale tenderebbe ad esplorare più velocemente.