Hamiltonowskie Monte Carlo

Hamiltona Monte Carlo (pierwotnie znany jako hybrydowy Monte Carlo ) to metoda Monte Carlo łańcucha Markowa służąca do uzyskiwania sekwencji losowych próbek , które zbiegają się w celu rozłożenia zgodnie z docelowym rozkładem prawdopodobieństwa, dla którego bezpośrednie próbkowanie jest trudne. Tej sekwencji można użyć do oszacowania całek w odniesieniu do docelowego rozkładu ( wartości oczekiwane ).

Hamiltonian Monte Carlo odpowiada przykładowi algorytmu Metropolisa-Hastingsa z ewolucją dynamiki Hamiltona symulowaną przy użyciu odwracalnego w czasie i zachowującego objętość integratora numerycznego (zwykle integratora przeskakującego ), aby zaproponować przejście do nowego punktu w przestrzeń stanów. W porównaniu z użyciem błądzenia losowego Gaussa rozkład propozycji w algorytmie Metropolisa-Hastingsa, hamiltonian Monte Carlo zmniejsza korelację między kolejnymi próbkowanymi stanami, proponując ruchy do odległych stanów, które zachowują wysokie prawdopodobieństwo akceptacji ze względu na przybliżone właściwości oszczędzania energii symulowanej dynamiki hamiltonowskiej przy użyciu integratora symplektycznego . Zmniejszona korelacja oznacza, że Monte Carlo potrzeba mniej próbek łańcucha Markowa .

Algorytm został pierwotnie zaproponowany przez Simona Duane'a, Anthony'ego Kennedy'ego, Briana Pendletona i Duncana Rowetha w 1987 roku do obliczeń w chromodynamice kwantowej kraty . W 1996 roku Radford M. Neal pokazał, w jaki sposób można zastosować tę metodę do szerszej klasy problemów statystycznych, w szczególności do sztucznych sieci neuronowych . Jednak ciężar konieczności dostarczania gradientów odpowiednich gęstości opóźnił szersze przyjęcie metody w statystyce i innych dyscyplinach ilościowych, aż w połowie 2010 roku twórcy Stan zaimplementowano HMC w połączeniu z automatycznym różnicowaniem .

Algorytm

Załóżmy, że docelowy rozkład na próbkę to ${\ Displaystyle f (\ mathbf {x})}$ dla ${\ Displaystyle \ mathbf {x} \ in \ mathbb {R} ^ {d}}$ ( ${\ Displaystyle d \ geq 1}$ ) i łańcuch próbek ${\ Displaystyle \ mathbf {X} _ {0}, \ mathbf {X} _ {1}, \ mathbf {X} _{2},\lddots }$ jest wymagane.

Równania Hamiltona są

${\ Displaystyle {\ Frac {{\ tekst {d}} x_ {i}} {{\ tekst {d}} t}} = {\ Frac {\ częściowe H} {\częściowe p_{i}}}}$

I

${\ Displaystyle {\ dfrac {{\ tekst {d}} p_ {i}} {{\ tekst {d}} t}} = - {\ dfrac {\ częściowe H}{\częściowe x_{i}}}}$

gdzie $,$ są odpowiednio $displaystyle$ składową wektora położenia $jest$ a p \ $}$ hamiltonianem . Niech będzie $macierzą$ mas , która jest symetryczna i określona, ​​to jest hamiltonian

${\ Displaystyle H (\ mathbf {x}, \ mathbf {p}) = U (\ mathbf {x}) + {\ dfrac {1}{2}}\mathbf {p} ^{\text{T}}M^{-1}\mathbf {p} }$

gdzie ${\ Displaystyle U (\ mathbf {x})}$ jest energią potencjalną . Energia potencjalna celu jest podana jako

${\ Displaystyle U (\ mathbf {x}) = - \ ln f (\ mathbf {x})}$

co wynika z czynnika Boltzmanna .

Algorytm wymaga dodatniej liczby całkowitej dla liczby kroków żaby przeskakującej $i$ dodatniej dla rozmiaru kroku ${\ displaystyle \ Delta t}$ . Załóżmy, że łańcuch jest w ${\ Displaystyle \ mathbf {X} _ {n} = \ mathbf {x} _ {n}}$ . Niech ${\ Displaystyle \ mathbf {x} _ {n} (0) = \ mathbf {x} _ {n}}$ . Najpierw losowy pęd Gaussa ${\ Displaystyle \ mathbf {p} _ {n} (0)}$ jest rysowany z ${\ Displaystyle {\ tekst {N}} \ lewo (\ mathbf {0}, M \ prawej)}$ . Następnie cząstka będzie działać przez czas pod dynamiką hamiltonowską $się$ to poprzez numeryczne rozwiązanie równań Hamiltona przy użyciu algorytmu przeskakującej Wektory położenia i pędu po czasie przy użyciu algorytmu żaby przeskakującej to ${\ displaystyle \ Delta t}$

${\ Displaystyle \ mathbf {p} _ {n} \ lewo (t + {\ dfrac {\ delta t} {2}} \ prawej) = \ mathbf {p} _ {n} (t )-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t) }}$

${\ Displaystyle \ mathbf {x} _ {n} (t + \ Delta t) = \ mathbf {x} _ {n }(t)+\Delta tM^{-1}\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}$

${\ Displaystyle \ mathbf {p} _ {n} (t + \ delta t) = \ mathbf {p} _ {n} \ lewo (t + {\ dfrac {\ delta t} {2}}\right)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _ {n}(t+\Delta t)}}$

Równania te mają być stosowane do ${\ Displaystyle \ mathbf {x} _ {n} (0)}$ i ${\ Displaystyle \ mathbf {p} _ {n} (0)}$${ \ Displaystyle L}$ razy, aby uzyskać ${\ Displaystyle \ mathbf {x} _ {n} (L \ Delta t)}$ i ${\ Displaystyle \ mathbf {p} _ {n}(L\Delta t)}$ .

Algorytm żaby przeskakującej jest przybliżonym rozwiązaniem ruchu nieoddziałujących cząstek klasycznych. Jeśli jest dokładne, rozwiązanie nigdy nie zmieni początkowego losowo wygenerowanego rozkładu energii, ponieważ energia jest zachowana dla każdej cząstki w obecności klasycznego pola energii potencjalnej. Aby osiągnąć rozkład równowagi termodynamicznej, cząstki muszą oddziaływać na przykład z otaczającą kąpielą cieplną, tak aby cały układ mógł przyjmować różne energie z prawdopodobieństwem zgodnym z rozkładem Boltzmanna.

Jednym ze sposobów przesunięcia układu w kierunku rozkładu równowagi termodynamicznej jest zmiana stanu cząstek za pomocą Metropolisa-Hastingsa . Więc najpierw stosuje się krok żaby przeskakującej, a następnie krok Metropolisa-Hastingsa.

Przejście od ${\ Displaystyle \ mathbf {X} _ {n} = \ mathbf {x} _ {n}}$ do ${\ Displaystyle \ mathbf {X} _ {n + 1} }$ jest

${\ Displaystyle \ mathbf {X} _ {n + 1} | \ mathbf {X} _ {n} = \ mathbf {x} _ {n} = {\ rozpocząć {przypadki} \ mathbf {x} _ {n} (L\Delta t)&{\text{z prawdopodobieństwem }}\alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right )\\\mathbf {x} _{n}(0)&{\text{inaczej}}\end{przypadki}}}$

Gdzie

$\ Displaystyle \ alfa \ lewo (\ mathbf {x} _ {n} (0), \ mathbf {x} _ {n} (L \ Delta t) \ prawej) = {\ text{min}}\left(1,{\dfrac {\exp \left[-H(\mathbf {x} _{n}(L\Delta t),\mathbf {p} _{n}(L\ Delta t))\right]}{\exp \left[-H(\mathbf {x} _{n}(0),\mathbf {p} _{n}(0))\right]}}\right )}$

Powtarza się to, aby uzyskać ${\ Displaystyle \ mathbf {X} _ {n + 1}, \ mathbf {X} _ {n + 2}, \ mathbf {X} _{n+3},\ldkropki}$ .

Brak samplera do zawracania

Próbnik bez zawracania (NUTS) jest rozszerzeniem polegającym na $automatycznym$ . Strojenie $krytyczne$ . Na przykład w przypadku jednowymiarowym potencjał wynosi ${\ Displaystyle U (x) = kx ^ {2}/2},$$\ sqrt {k}})}$ co odpowiada potencjałowi prostego oscylatora harmonicznego . Dla $zbyt dużych cząstek$ będzie oscylować, a tym samym tracić czas obliczeniowy. Dla $małej$ cząstki cząstka będzie zachowywać się jak błądzenie losowe.

Luźno, NUTS losowo uruchamia dynamikę Hamiltona zarówno do przodu, jak i do tyłu w czasie, aż spełniony jest warunek zawracania. Kiedy tak się dzieje, dla próbki MCMC wybierany jest losowy punkt ze ścieżki i proces jest powtarzany od tego nowego punktu.

Mówiąc szczegółowo, drzewo binarne jest konstruowane w celu śledzenia ścieżki kroków żaby przeskakującej. Aby wytworzyć próbkę MCMC, przeprowadza się procedurę iteracyjną. ${\ Displaystyle U_ {n} \ sim {\ tekst {jednolity}} (0, \ exp (-$ Jednolity jest próbkowane. Niech ${\ Displaystyle \ mathbf {x} _ {n}$$^ {+}}$ odpowiednio i pędem cząstki przedniej. Podobnie, ${\ Displaystyle \ mathbf {x} _ {n} ^ {-}}$ i ${\ Displaystyle \ mathbf {p} _ {n} ^ {-}}$ dla wstecznej cząstki. W każdej iteracji drzewo binarne wybiera losowo i równomiernie, aby przesunąć cząstkę do przodu w czasie lub cząstkę do tyłu w czasie. Również dla każdej iteracji liczba kroków żaby przeskakującej wzrasta o współczynnik 2. Na przykład w pierwszej iteracji cząstka do przodu porusza się do przodu w czasie za pomocą 1 kroku żaby przeskakującej. W następnej iteracji cząsteczka wsteczna cofa się w czasie, wykonując 2 kroki żaby przeskakującej.

Procedura iteracyjna jest kontynuowana do momentu spełnienia warunku zawracania, tj

${\ Displaystyle (\ mathbf {x} _ {n} ^ {+} - \ mathbf { x} _{n}^{-})\cdot \mathbf {p} _{n}^{-}<0\quad {\text{or}}\quad (\mathbf {x} _{n}^ {+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{+}<0}$

lub gdy hamiltonian staje się niedokładny

${\ Displaystyle \ exp \ lewo [-H (\ mathbf {x} _ {n} ^ {+}, \ mathbf {p} _{n}^{+})+\delta \right]<U_{n}}$

Lub

${\ Displaystyle \ exp \ lewo [-H (\ mathbf {x} _ {n} ^ {-}, \ mathbf {p} _{n}^{-})+\delta \right]<U_{n}}$

gdzie na przykład ${\ displaystyle \ delta = 1000}$ .

Po spełnieniu warunku zawracania następna próbka MCMC , $jest$ przeskakującej wytyczonej przez drzewo ${\ Displaystyle \ {\ mathbf {x} _ {n} ^ {-}, \ ldots, \ mathbf {x} _ {n} (- \ Delta t), \ mathbf {x} _ {n} (0) ,\mathbf {x} _{n}(\Delta t),\ldots ,\mathbf {x} _{n}^{+}\}} który$ spełnia

${\ Displaystyle U_ {n}< \ exp \ lewo [-H (\ mathbf {x_ {n + 1}} \ mathbf { p_{n+1})} \prawo]}$

Zwykle jest to spełnione, jeśli pozostałe parametry konsoli HMC są rozsądne.

Zobacz też

• Dynamiczna metoda Monte Carlo
• Oprogramowanie do modelowania molekularnego metodą Monte Carlo
• Stan
• Algorytm Langevina dostosowany do Metropolisa

Dalsza lektura

•   Neal, Radford M (2011). „MCMC przy użyciu dynamiki hamiltonowskiej” (PDF) . W Steve'ie Brooksie; Andrzeja Gelmana; Galin L. Jones; Xiao-Li Meng (red.). Podręcznik łańcucha Markowa Monte Carlo . Chapmana i Halla/CRC. ISBN 9781420079418 .
• Betancourt, Michael (2018). „Koncepcyjne wprowadzenie do hamiltonowskiego Monte Carlo”. ar Xiv : 1701.02434 .
•   Barbu, Adrian; Zhu, Song-Chun (2020). „Hamiltonian i Langevin Monte Carlo”. Metody Monte Carlo . Singapur: Springer. s. 281–326. ISBN 978-981-13-2970-8 .

Linki zewnętrzne

• Betancourt, Michael. „Efektywne wnioskowanie bayesowskie z hamiltonowskim Monte Carlo” . MLSS Islandia 2014 – przez YouTube .
• McElreath, Richard. „Łańcuch Markowa Monte Carlo” . Statystyczne przemyślenie 2022 – za pośrednictwem YouTube .
• Hamiltonowskie Monte Carlo od podstaw
• Metody optymalizacji i Monte Carlo