Algorytm Metropolisa-Hastingsa
Algorytm Metropolisa-Hastingsa to technika statystyczna klasy MCMC (próbkowanie Monte Carlo łańcuchami Markowa), umożliwiająca stochastyczne oszacowanie całek (takich jak estymatory) oraz rozkładów prawdopodobieństwa dla skomplikowanych systemów, które są zbyt złożone do analizy analitycznej, jak na przykład układy wielowymiarowe. To narzędzie znajduje zastosowanie między innymi w wnioskowaniu bayesowskim oraz modelowaniu systemów fizycznych. Po raz pierwszy procedura ta została opisana publicznie przez zespół Metropolisa w 1953 roku, a jej rozwinięcie miało miejsce w 1970 roku dzięki Hastingsowi. Metody Monte Carlo tego typu były również wykorzystywane przez Enrico Fermiego i Stanisława Ulama podczas tajnych prac nad Projektem Manhattan.
Systemy o wielu wymiarach są obarczone zjawiskiem przekleństwa wymiarowości, co oznacza, że wraz ze wzrostem liczby wymiarów problemu, liczba wymaganych obserwacji do oszacowania ich właściwości rośnie w sposób wykładniczy. Metody Monte Carlo oparte na łańcuchach Markowa są w dużym stopniu odporne na ten problem, ponieważ nie wymagają analitycznego rozwiązania ani przeszukiwania całej przestrzeni problemu. Zamiast tego, stosują iteracyjne próbkowanie stochastyczne, które z każdą iteracją skupia się coraz bardziej na centralnych obszarach rozkładów.
Wprowadzenie metod MCMC oraz dostępność komputerów umożliwiły przezwyciężenie wielu ograniczeń związanych z analizą złożonych problemów w naukach empirycznych, które wcześniej opierały się głównie na uproszczonych metodach tradycyjnej statystyki. Dokładność uzyskanych oszacowań w nowych metodach zależy od dopasowania parametrów łańcucha Markowa do oczekiwanego rozkładu oraz od liczby powtórzeń próbkowania.
Algorytm Metropolisa-Hastingsa przedstawia jeden z prostych w implementacji sposobów na skonstruowanie łańcucha Markowa za pomocą błądzenia losowego, co pozwala na rozwiązanie wielu typowych problemów. W późniejszych latach opracowano również metody dostosowane do szerszych klas problemów, na przykład śledzące kształt rozkładu z wykorzystaniem wektorów pędu, szczególnie hamiltonianów.
Algorytm
Algorytm opiera się na powtarzaniu coraz dokładniejszych iteracji próbkowania rozkładu P(x) w oparciu o błądzenie losowe według wybranej funkcji g, zgodnie z poniżej przedstawionymi krokami. Wymaga on, aby próbkowany rozkład był stabilny i nieokresowy, ponieważ w przeciwnym razie łańcuch Markowa może ulec zapętleniu.
1. Losowa inicjalizacja
Pierwszy stan x łańcucha Markowa jest wybierany losowo z obszaru danego rozkładu losowego.
2. Wybranie następnego stanu x′
Następny stan jest wybierany zgodnie z funkcją błądzenia losowego g(x’|x), określoną przez użytkownika. Funkcja g może być na przykład oparta na rozkładzie normalnym lub rozkładzie t Studenta (a w przypadku funkcji wielu zmiennych, na wielowymiarowym rozkładzie normalnym lub wielowymiarowym rozkładzie t Studenta).
3. Przyjęcie lub odrzucenie stanu x′
Ocena nowego stanu odbywa się według kryterium wybranego przez użytkownika. Na przykład, Metropolis proponuje funkcję:
A(x’|x) = min(1, (P(x’) / P(x)) * (g(x|x’) / g(x’|x))). Akceptuje się przejście do nowego stanu x’, jeśli proponowane parametry pozwalają na oszacowanie rozkładu z wyższym prawdopodobieństwem.
4. Powrót do punktu 2
Testowanie lokalnych stanów x’ odbywa się przy wybranej przez użytkownika liczbie powtórzeń T.
5. Przyjęcie nowego stanu
Po przetestowaniu stanu x’, jest on zapisywany jako nowy x, a algorytm poszukuje kolejnego punktu x’.
Przykład próbkowania
Poniżej znajduje się kod w języku Python, który ilustruje próbkowanie rozkładu normalnego, nieunormowanego do 1, przy użyciu metody Metropolisa. Dodatkowo, zawiera wykresy porównujące histogram uzyskany z próbkowania z unormowanym rozkładem normalnym. Program można uruchomić, korzystając na przykład z darmowego notatnika Google Colab online.
Oprogramowanie dla metod MCMC
Popularną implementacją algorytmu Metropolisa-Hastingsa oraz innych zaawansowanych metod MCMC jest darmowe i otwarte oprogramowanie STAN, które jest dostępne w otwartym pakiecie statystycznym R i wykorzystywane w metodach wnioskowania bayesowskiego.
