Bagaimana cara kerja Hamiltonian Monte Carlo?

Question

Saya membuat grafik di bawah ini untuk menjelaskan bagaimana saya saat ini memahami algoritma HMC. Saya ingin verifikasi dari ahli materi pelajaran jika pemahaman ini benar atau tidak. Teks di slide di bawah ini disalin di bawah untuk kemudahan akses:

---

Hamiltonian Monte Carlo: Sebuah satelit mengorbit planet. Semakin dekat satelit ke planet, semakin besar efek gravitasi. Ini berarti, (A) energi potensial yang lebih tinggi dan (B) energi kinetik yang lebih tinggi diperlukan untuk mempertahankan orbit. Energi kinetik yang sama pada jarak yang lebih jauh dari planet ini akan mengeluarkan satelit dari orbit. Satelit tersebut bertugas mengumpulkan foto dari kawasan geografis tertentu. Semakin dekat satelit mengorbit planet, semakin cepat ia bergerak di orbit, semakin sering ia melewati wilayah tersebut, semakin banyak foto yang dikumpulkannya. Sebaliknya, semakin jauh sebuah satelit dari planet ini, semakin lambat ia bergerak di orbit, semakin sedikit ia melewati wilayah tersebut, semakin sedikit foto yang dikumpulkannya. Dalam konteks pengambilan sampel, jarak dari planet mewakili jarak dari ekspektasi distribusi. Area dengan kemungkinan rendah jauh dari harapan; ketika “mengorbit kemungkinan ini,” energi kinetik yang lebih rendah berarti lebih sedikit sampel yang dikumpulkan selama interval waktu yang tetap, sedangkan ketika mengorbit kemungkinan yang lebih tinggi berarti lebih banyak sampel yang dikumpulkan dengan interval waktu tetap yang sama. Dalam orbit tertentu, energi total, kinetik dan potensial, konstan; Namun, hubungan keduanya tidak sederhana. Persamaan Hamiltonian menghubungkan perubahan satu dengan lainnya. Yakni, gradien posisi terhadap waktu sama dengan momentum. Dan gradien momentum terhadap waktu sama dengan gradien energi potensial terhadap posisi. Untuk menghitung sejauh mana satelit akan menempuh jalur orbitnya, integrasi lompatan harus digunakan, yang secara berulang memperbarui momentum dan vektor posisi. Dalam konteks pengambilan sampel, kemungkinan dianalogikan dengan jarak dari planet dan gradien energi potensial sehubungan dengan posisi adalah gradien dari fungsi kepadatan probabilitas sehubungan dengan parameter masukannya, x. Informasi ini memungkinkan jalur orbit di sekitar berbagai masukan, X, yang sesuai dengan kemungkinan yang sama, y, untuk dieksplorasi.
Namun, kami tidak hanya tertarik untuk menjelajahi satu kemungkinan, kami harus menjelajahi beberapa jalur orbit. Untuk mencapai ini, momentum harus diperbesar secara acak, membawa satelit lebih dekat atau lebih jauh dari planet ini. "Tendangan momentum" acak ini memungkinkan kemungkinan orbit yang berbeda. Untungnya, persamaan hamiltonian memastikan bahwa tidak peduli kemungkinannya, jumlah sampel yang dikumpulkan sebanding dengan kemungkinannya, sehingga sampel yang dikumpulkan mengikuti bentuk distribusi target.

---

Pertanyaan saya adalah - Apakah ini cara yang akurat untuk berpikir tentang cara kerja Hamiltonian Monte Carlo?

Edit:

Saya telah menerapkan dalam beberapa kode berdasarkan pemahaman saya tentang algoritme. Ia bekerja untuk gaussian dengan mu = 0, sigma = 1. Tetapi jika saya mengubah sigma itu rusak. Semua wawasan akan dihargai.

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

Sekarang berfungsi di sini:

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()

Tapi itu rusak ketika saya mengubah sigma menjadi 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()

Ada ide? Saya merasa seperti saya hampir "mendapatkannya".

Answer 1

Sebelum menjawab pertanyaan tentang cara intuitif untuk berpikir tentang Hamiltonian Monte Carlo, mungkin yang terbaik adalah memahami MCMC biasa dengan benar. Mari kita kesampingkan metafora satelit untuk saat ini.

MCMC berguna ketika Anda menginginkan sampel unbiassed dari distribusi di mana Anda hanya memiliki sesuatu yang tersedia yang sebanding dengan PDF, tetapi bukan PDF itu sendiri. Ini muncul dalam (misalnya) simulasi fisika: PDF diberikan oleh distribusi Boltzmann, p ~ exp (-E / kT), tetapi hal yang dapat Anda hitung untuk konfigurasi sistem apa pun adalah E, bukan p. Konstanta proporsionalitas tidak diketahui, karena integral dari exp (-E / kT) pada seluruh ruang dari konfigurasi yang memungkinkan biasanya terlalu sulit untuk dihitung. MCMC memecahkan masalah tersebut dengan melakukan jalan acak dengan cara tertentu, di mana probabilitas untuk mengambil ("menerima") setiap langkah terkait dengan rasio nilai p (konstanta proporsionalitas dibatalkan). Seiring waktu, distribusi sampel yang diterima dari random walk menyatu dengan PDF yang kita inginkan, tanpa perlu menghitung p secara eksplisit.

Perhatikan bahwa di atas, metode apa pun untuk mengambil langkah acak sama validnya, selama pejalan acak dapat menjelajahi seluruh ruang. Kriteria penerimaan menjamin bahwa sampel yang dipilih menyatu dengan PDF asli. Dalam praktiknya, distribusi gaussian di sekitar sampel saat ini digunakan (dan sigma dapat divariasikan sehingga fraksi langkah yang diterima tetap relatif tinggi). Pada prinsipnya, tidak ada yang salah dengan mengambil langkah-langkah dari distribusi kontinu lainnya ("distribusi lompatan") di sekitar sampel saat ini, meskipun konvergensi mungkin jauh lebih lambat.

Sekarang, Hamiltonian Monte Carlo memperluas metafora fisika dengan secara khusus mencoba mengambil langkah ke arah yang lebih mungkin diterima daripada langkah gaussian. Langkah-langkah yang akan diambil oleh integrator lompatan, jika mencoba menyelesaikan gerakan sistem yang energi potensinya E. Persamaan gerak ini juga mencakup istilah energi kinetik, dengan "massa" (bukan fisik) dan "momentum". Langkah-langkah yang dilakukan integrator lompatan dalam "waktu" kemudian diteruskan sebagai proposal ke algoritme MCMC.

Mengapa ini berhasil? MC gaussian mengambil langkah-langkah dengan jarak yang sama ke segala arah dengan probabilitas yang sama; satu-satunya hal yang membiaskannya ke area PDF yang lebih padat penduduknya adalah bahwa langkah-langkah ke arah yang salah kemungkinan besar akan ditolak. Hamiltonian MC mengusulkan langkah-langkah baik dalam arah gradien E, dan arah gerakan yang terakumulasi dalam langkah-langkah terakhir (arah dan besaran "momentum"). Hal ini memungkinkan eksplorasi ruang yang lebih cepat, dan juga kemungkinan yang lebih tinggi untuk menjangkau wilayah yang lebih padat dengan lebih cepat.

Sekarang, metafora satelit: Saya pikir ini bukan cara yang berguna untuk memikirkannya. Satelit bergerak tepat pada orbitnya; apa yang Anda miliki di sini cukup acak, lebih seperti partikel gas dalam wadah dengan partikel lain. Setiap tabrakan acak memberi Anda "langkah"; seiring waktu partikel akan berada di mana-mana dalam wadah dengan probabilitas yang sama (karena PDF di sini sama di mana-mana, kecuali dinding yang mewakili energi yang sangat tinggi / efektif nol PDF). Gaussian MCMC seperti partikel bermassa nol efektif yang berjalan secara acak (atau partikel bermassa bukan-nol dalam medium yang relatif kental): ia akan sampai ke sana melalui gerakan brownian, tetapi tidak harus cepat. Hamiltonian MC adalah partikel bermassa bukan nol: ia dapat mengumpulkan momentum yang cukup untuk terus bergerak ke arah yang sama meskipun terjadi tabrakan, dan karenanya kadang-kadang dapat menembak dari satu ujung wadah ke ujung lain (tergantung pada massa vs frekuensi / besarnya tabrakan). Itu masih akan memantul dari dinding, tentu saja, tetapi secara umum cenderung menjelajah lebih cepat.