Algorytm Metropolis-Hastings - Metropolis–Hastings algorithm

Propozycja dystrybucja Q proponuje następny punkt, do którego błądzenia losowego może przenieść.

W statystykach i fizyce statystycznych The Metropolisa Hastings algorytm jest Monte Carlo łańcuchów Markowa Sposób (MCMC) dla uzyskania sekwencji losowych próbek z rozkładu prawdopodobieństwa , z których bezpośrednie pobieranie jest trudne. Ta sekwencja może byćużyta do przybliżenia rozkładu (np. do wygenerowania histogramu ) lub do obliczenia całki (np.wartości oczekiwanej). Metropolis-Hastings i inne algorytmy MCMC są zwykle używane do próbkowania z rozkładów wielowymiarowych, zwłaszcza gdy liczba wymiarów jest duża. W przypadku rozkładów jednowymiarowych istnieją zwykle inne metody (np. próbkowanie adaptacyjne z odrzuceniem ), które mogą bezpośrednio zwracać niezależne próbki z rozkładu i są one wolne od problemu autokorelacji próbek, który jest nieodłączny w metodach MCMC.

Historia

Algorytm został nazwany na cześć Nicholasa Metropolisa , który wraz z Arianną W. Rosenbluth , Marshallem Rosenbluthem , Augustą H. Teller i Edwardem Tellerem napisał w 1953 r. artykuł Equation of State Calculations by Fast Computing Machines . W tym artykule zaproponowano algorytm dla przypadku symetrycznych rozkładów propozycji, a WK Hastings rozszerzył go na bardziej ogólny przypadek w 1970 roku.

Pewne kontrowersje budzi zasługa opracowania algorytmu. Metropolis, obeznany z obliczeniowymi aspektami metody, ukuł we wcześniejszym artykule ze Stanisławem Ulamem termin „Monte Carlo” i kierował grupą w Zakładzie Teoretycznym, która zaprojektowała i zbudowała komputer MANIAC I używany w eksperymentach w 1952. Jednak przed 2003 r. nie było szczegółowego opisu rozwoju algorytmu. Następnie, na krótko przed śmiercią, Marshall Rosenbluth wziął udział w konferencji w LANL w 2003 roku z okazji 50. rocznicy publikacji z 1953 roku. Na tej konferencji Rosenbluth opisał algorytm i jego rozwój w prezentacji zatytułowanej „Genesis of the Monte Carlo Algorithm for Statistical Mechanics”. Dalsze wyjaśnienia historyczne zostały zawarte przez Gubernatisa w artykule z 2005 roku, opowiadającym o 50-leciu konferencji. Rosenbluth wyjaśnia, że ​​on i jego żona Arianna wykonali tę pracę, a Metropolis nie odegrało żadnej roli w rozwoju poza zapewnieniem czasu komputerowego.

Jest to sprzeczne z relacją Edwarda Tellera, który w swoich pamiętnikach stwierdza, że ​​pięciu autorów artykułu z 1953 roku pracowało razem przez „dni (i noce)”. W przeciwieństwie do tego, szczegółowa relacja Rosenblutha przypisuje Tellerowi kluczową, ale wczesną sugestię, aby „wykorzystać mechanikę statystyczną i wziąć średnie zespołowe zamiast podążać za szczegółową kinematykę”. To, mówi Rosenbluth, skłoniło go do myślenia o uogólnionym podejściu Monte Carlo – temacie, który, jak twierdzi, często omawiał z Johnem Von Neumannem . Arianna Rosenbluth opowiedziała (w Gubernatis w 2003 r.), że Augusta Teller rozpoczęła pracę z komputerem, ale sama Arianna przejęła go i napisała kod od zera. W ustnej historii zarejestrowanej na krótko przed śmiercią Rosenbluth ponownie przypisuje Tellerowi postawienie pierwotnego problemu, siebie za rozwiązanie go, a Ariannie za programowanie komputera. Jeśli chodzi o reputację, nie ma powodu, aby kwestionować konto Rosenblutha. W biografii Rosenblutha Freeman Dyson pisze:

Wielokrotnie przychodziłem do Rosenbluth, zadając mu pytanie [...] i otrzymując odpowiedź w ciągu dwóch minut. Wtedy zwykle potrzebowałem tygodnia ciężkiej pracy, aby zrozumieć szczegółowo, dlaczego odpowiedź Rosenblutha była słuszna. Miał niesamowitą zdolność przejrzenia skomplikowanej sytuacji fizycznej i znalezienia właściwej odpowiedzi za pomocą fizycznych argumentów. Enrico Fermi był jedynym znanym mi fizykiem, który dorównywał Rosenbluthowi w swoim intuicyjnym pojmowaniu fizyki.

Intuicja

Algorytm Metropolis-Hastings może pobierają próbki z dowolnego rozkładu prawdopodobieństwa z gęstością prawdopodobieństwa , pod warunkiem, że wiemy Funkcja A proporcjonalny do gęstości i wartości mogą być obliczone. Wymóg, który musi być tylko proporcjonalny do gęstości, a nie dokładnie równy, czyni algorytm Metropolisa-Hastingsa szczególnie użytecznym, ponieważ obliczenie niezbędnego współczynnika normalizacji jest często niezwykle trudne w praktyce.

Algorytm Metropolis-Hastings działa, generując sekwencję wartości próbek w taki sposób, że w miarę wytwarzania coraz większej liczby wartości próbek rozkład wartości jest bardziej zbliżony do pożądanego rozkładu. Te wartości próbek są tworzone iteracyjnie, przy czym rozkład następnej próbki jest zależny tylko od bieżącej wartości próbki (co powoduje, że sekwencja próbek tworzy łańcuch Markowa ). W szczególności w każdej iteracji algorytm wybiera kandydata na następną wartość próbki na podstawie bieżącej wartości próbki. Następnie, z pewnym prawdopodobieństwem, kandydat jest albo akceptowany (wówczas wartość kandydująca jest używana w następnej iteracji) lub odrzucany (wówczas kandydująca wartość jest odrzucana, a aktualna wartość jest ponownie używana w następnej iteracji) — prawdopodobieństwo akceptacji jest określana przez porównanie wartości funkcji bieżącej i kandydującej wartości próbki w odniesieniu do pożądanego rozkładu.

W celu ilustracji poniżej opisano algorytm Metropolis, szczególny przypadek algorytmu Metropolis-Hastings, w którym funkcja propozycji jest symetryczna.

Algorytm Metropolis (symetryczny rozkład propozycji)

Niech będzie funkcją proporcjonalną do pożądanej funkcji gęstości prawdopodobieństwa (inaczej rozkład docelowy).

  1. Inicjalizacja: Wybierz dowolny punkt, który będzie pierwszą obserwacją w próbce i wybierz dowolną gęstość prawdopodobieństwa (czasami zapisywaną ), która sugeruje kandydata do następnej wartości próbki , biorąc pod uwagę poprzednią wartość próbki . W tej sekcji zakłada się, że jest symetryczny; innymi słowy, musi spełniać . Zwykły wybór to niech będzie rozkładem Gaussa skoncentrowane na , tak, że punkty bliżej są bardziej skłonni do odwiedzin obok, zrobienie sekwencji próbek do błądzenia losowego . Funkcja ta nazywana jest gęstością propozycji lub dystrybucją skokową .
  2. Dla każdej iteracji t :
    • Wygeneruj kandydata do następnej próbki, wybierając z dystrybucji .
    • Oblicz ten współczynnik akceptacji , który będzie używany do zdecydować, czy przyjąć lub odrzucić kandydata. Ponieważ f jest proporcjonalne do gęstości P , mamy to .
    • Zaakceptuj lub odrzuć :
      • Wygeneruj jednolitą liczbę losową .
      • Jeśli , to zaakceptuj kandydata poprzez ustawienie ,
      • Jeśli , to odrzuć kandydata i ustaw zamiast tego.

Algorytm ten działa poprzez losowe próby poruszania się po przestrzeni próbki, czasami akceptując ruchy, a czasami pozostając na miejscu. Należy zauważyć, że współczynnik akceptacji wskazuje, jak prawdopodobne jest nowo proponowana próbka w odniesieniu do próbki bieżącej, zgodnie z rozkładem, którego gęstość wynosi . Jeśli spróbujemy przenieść się do punktu, który jest bardziej prawdopodobny niż istniejący punkt (tj. punkt w regionie o większej gęstości odpowiadającego ), zawsze zaakceptujemy ruch. Jeśli jednak spróbujemy przejść do mniej prawdopodobnego punktu, czasami odrzucimy ruch, a im większy względny spadek prawdopodobieństwa, tym większe prawdopodobieństwo, że odrzucimy nowy punkt. W związku z tym będziemy mieć tendencję do pozostawania w (i zwracania dużej liczby próbek z) regionów o wysokiej gęstości , podczas gdy tylko sporadycznie odwiedzamy regiony o niskiej gęstości. Intuicyjnie dlatego ten algorytm działa i zwraca próbki, które podążają za pożądanym rozkładem z gęstością .

W porównaniu z algorytmem takim jak adaptacyjne próbkowanie odrzucania, który bezpośrednio generuje niezależne próbki z rozkładu, Metropolis-Hastings i inne algorytmy MCMC mają szereg wad:

  • Próbki są skorelowane. Nawet jeśli na dłuższą metę prawidłowo podążają , zestaw pobliskich próbek będzie ze sobą skorelowany i nie będzie poprawnie odzwierciedlał rozkładu. Oznacza to, że jeśli chcemy mieć zestaw niezależnych próbek, musimy wyrzucić większość próbek i pobierać tylko każdą n- tą próbkę, dla pewnej wartości n (zwykle określanej przez badanie autokorelacji między sąsiednimi próbkami). Autokorelację można zmniejszyć, zwiększając szerokość skoku (średnią wielkość skoku, która jest związana z wariancją rozkładu skoków), ale zwiększy to również prawdopodobieństwo odrzucenia proponowanego skoku. Zbyt duży lub zbyt mały skokowy rozmiar doprowadzi do wolno mieszającego się łańcucha Markowa, tj. wysoce skorelowanego zestawu próbek, tak że do uzyskania rozsądnego oszacowania dowolnej pożądanej właściwości rozkładu potrzebna będzie bardzo duża liczba próbek.
  • Chociaż łańcuch Markowa ostatecznie zbiega się do pożądanego rozkładu, początkowe próbki mogą mieć bardzo różny rozkład, zwłaszcza jeśli punkt początkowy znajduje się w regionie o niskiej gęstości. W rezultacie, wypalenia zazwyczaj konieczne jest czas, w którym początkowa liczba próbek (np pierwszym 1000 lub więcej) są brane pod uwagę.

Z drugiej strony, większość prostych metod próbkowania odrzuceń cierpi z powodu „ przekleństwa wymiarowości ”, gdzie prawdopodobieństwo odrzucenia wzrasta wykładniczo w funkcji liczby wymiarów. Metropolis-Hastings, podobnie jak inne metody MCMC, nie mają tego problemu w takim stopniu, a zatem są często jedynymi dostępnymi rozwiązaniami, gdy liczba wymiarów rozkładu, które mają być próbkowane, jest duża. W rezultacie metody MCMC są często wybieranymi metodami wytwarzania próbek z hierarchicznych modeli bayesowskich i innych wysokowymiarowych modeli statystycznych stosowanych obecnie w wielu dyscyplinach.

W rozkładach wielowymiarowych klasyczny algorytm Metropolis-Hastings opisany powyżej obejmuje wybór nowego wielowymiarowego punktu próbkowania. Gdy liczba wymiarów jest duża, znalezienie odpowiedniego rozkładu przeskoków do użycia może być trudne, ponieważ poszczególne wymiary zachowują się w bardzo różny sposób, a szerokość przeskoku (patrz powyżej) musi być „właściwa” dla wszystkich wymiarów jednocześnie, aby unikać zbyt powolnego mieszania. Alternatywne podejście, które często działa lepiej w takich sytuacjach, znane jako próbkowanie Gibbsa , polega na wybraniu nowej próbki dla każdego wymiaru oddzielnie od pozostałych, zamiast wybierania próbki dla wszystkich wymiarów naraz. W ten sposób problem próbkowania z przestrzeni potencjalnie wielowymiarowej zostanie zredukowany do zbioru problemów do próbkowania z małej wymiarowości. Ma to szczególne zastosowanie, gdy rozkład wielowymiarowy składa się z zestawu pojedynczych zmiennych losowych, w których każda zmienna jest uwarunkowana jedynie niewielką liczbą innych zmiennych, jak ma to miejsce w większości typowych modeli hierarchicznych . Poszczególne zmienne są następnie próbkowane pojedynczo, przy czym każda zmienna jest uwarunkowana najnowszymi wartościami wszystkich pozostałych. Różne algorytmy mogą być użyte do wyboru tych pojedynczych próbek, w zależności od dokładnej postaci rozkładu wielowymiarowego: niektóre możliwości to metody próbkowania adaptacyjnego odrzucenia, algorytm próbkowania adaptacyjnego odrzucenia Metropolis, prosty jednowymiarowy krok Metropolis-Hastings lub próbkowanie plastra .

Formalne pochodzenie

Celem algorytmu Metropolis-Hastings jest wygenerowanie zbioru stanów zgodnie z pożądanym rozkładem . Aby to osiągnąć, algorytm wykorzystuje proces Markowa , który asymptotycznie osiąga unikalny rozkład stacjonarny, taki jak .

Proces Markowa jest jednoznacznie określony przez jego prawdopodobieństwa przejścia , prawdopodobieństwo przejścia z dowolnego danego stanu do dowolnego innego danego stanu . Posiada unikalny rozkład stacjonarny, gdy spełnione są dwa warunki:

  1. Istnienie dystrybucji stacjonarnej : musi istnieć dystrybucja stacjonarna . Wystarczającym, ale nie koniecznym warunkiem jest szczegółowa równowaga , która wymaga, aby każde przejście było odwracalne: dla każdej pary stanów prawdopodobieństwo bycia w stanie i przejścia do stanu musi być równe prawdopodobieństwu bycia w stanie i przejścia do stanu , .
  2. Wyjątkowość dystrybucji stacjonarnej : dystrybucja stacjonarna musi być niepowtarzalna. Gwarantuje to ergodyczność procesu Markowa, który wymaga, aby każdy stan był (1) aperiodyczny — system nie powraca do tego samego stanu w stałych odstępach czasu; oraz (2) być dodatnie powtarzalne – oczekiwana liczba kroków powrotu do tego samego stanu jest skończona.

Algorytm Metropolis-Hastings polega na zaprojektowaniu procesu Markowa (poprzez skonstruowanie prawdopodobieństw przejścia), który spełnia dwa powyższe warunki, tak że jego rozkład stacjonarny jest wybrany jako . Wyprowadzenie algorytmu rozpoczyna się od warunku bilansu szczegółowego :

który jest przepisany jako

Podejście polega na podzieleniu przejścia na dwa podetapy; propozycję i akceptację-odrzucenie. Rozkład propozycja jest prawdopodobieństwo warunkowe proponując stan daną oraz rozkład akceptacja jest prawdopodobieństwem przyjąć proponowaną stan . Prawdopodobieństwo przejścia można zapisać jako iloczyn ich:

Wstawiając tę ​​relację do poprzedniego równania, mamy

Następnym krokiem w wyprowadzeniu jest wybór współczynnika akceptacji, który spełnia powyższy warunek. Jednym z powszechnych wyborów jest wybór Metropolis:

Dla tego współczynnika akceptacji Metropolis , albo albo i tak, warunek jest spełniony.

Algorytm Metropolis-Hastings składa się zatem z następujących elementów:

  1. Zainicjuj
    1. Wybierz stan początkowy .
    2. Ustaw .
  2. Powtarzać
    1. Wygeneruj losowy stan kandydujący zgodnie z .
    2. Oblicz prawdopodobieństwo akceptacji .
    3. Zaakceptuj lub odrzuć :
      1. wygeneruj jednolitą liczbę losową ;
      2. jeśli , to zaakceptuj nowy stan i ustaw ;
      3. if , odrzuć nowy stan i skopiuj stary stan do przodu .
    4. Przyrost : zestaw .

O ile zostaną spełnione określone warunki, empiryczny rozkład zapisanych stanów zbliży się do . Liczba iteracji ( ) wymaganych do skutecznego oszacowania zależy od liczby czynników, w tym relacji między dystrybucją propozycji i pożądanej dokładności oszacowania. Rozkład na dyskretnych przestrzeniach stanów musi być rzędu czasu autokorelacji procesu Markowa.

Należy zauważyć, że w ogólnym zagadnieniu nie jest jasne, jakiego rozkładu należy użyć, ani liczby iteracji niezbędnych do prawidłowego oszacowania; oba są dowolnymi parametrami metody, które należy dostosować do konkretnego problemu.

Użyj w całkowaniu numerycznym

Powszechnym zastosowaniem algorytmu Metropolis-Hastings jest obliczenie całki. W szczególności rozważ spację i rozkład prawdopodobieństwa nad , . Metropolis-Hastings może oszacować całkę postaci

gdzie jest (mierzalną) funkcją zainteresowania.

Rozważmy na przykład statystykę i jej rozkład prawdopodobieństwa , który jest rozkładem marginalnym . Załóżmy, że celem jest oszacowanie dla na ogonie . Formalnie można zapisać jako

a zatem oszacowanie można przeprowadzić przez oszacowanie oczekiwanej wartości funkcji wskaźnika , która wynosi 1, gdy i zero w przeciwnym razie. Ponieważ znajduje się na ogonie , prawdopodobieństwo narysowania stanu z ogonem jest proporcjonalne do , które z definicji jest małe. Algorytm Metropolis-Hastings może być tutaj użyty do próbkowania (rzadkich) stanów z większym prawdopodobieństwem, a tym samym do zwiększenia liczby próbek używanych do oszacowania ogonów. Można to zrobić np. stosując rozkład próbkowania faworyzujący te stany (np. za pomocą ).

Instrukcje krok po kroku

Załóżmy, że najnowsza próbkowana wartość to . Aby postępować zgodnie z algorytmem Metropolis-Hastings, rysujemy nowy stan propozycji z gęstością prawdopodobieństwa i obliczamy wartość

gdzie

jest prawdopodobieństwem (np. bayesowskim a posteriori) między proponowaną próbą a poprzednią próbą , a

to stosunek gęstości propozycji w dwóch kierunkach (od do i odwrotnie). Jest to równe 1, jeśli gęstość propozycji jest symetryczna. Następnie nowy stan jest wybierany zgodnie z następującymi zasadami.

Gdyby
w przeciwnym razie:

Łańcuch Markowa rozpoczyna się od dowolnej wartości początkowej , a algorytm jest uruchamiany przez wiele iteracji, aż ten stan początkowy zostanie „zapomniany”. Te próbki, które są wyrzucane, nazywane są wypaleniem . Pozostały zestaw przyjętych wartości reprezentują próbkę z rozkładu .

Algorytm działa najlepiej, gdy gęstość propozycji odpowiada kształtowi rozkładu docelowego , z którego bezpośrednie próbkowanie jest utrudnione, czyli . Jeśli używana jest gęstość propozycji Gaussa , parametr wariancji musi być dostrojony w okresie wygrzewania. Zwykle odbywa się to poprzez obliczenie współczynnika akceptacji , który jest ułamkiem proponowanych próbek, który jest akceptowany w oknie ostatnich próbek. Pożądany współczynnik akceptacji zależy od docelowego rozkładu, jednak teoretycznie wykazano, że idealny współczynnik akceptacji dla jednowymiarowego rozkładu Gaussa wynosi około 50%, spadając do około 23% dla dwuwymiarowego rozkładu Gaussa.

Jeśli jest za mały, łańcuch będzie się powoli mieszać (tj. wskaźnik akceptacji będzie wysoki, ale kolejne próbki będą powoli poruszać się po przestrzeni, a łańcuch będzie się powoli zbiegał do ). Z drugiej strony, jeśli jest zbyt duży, wskaźnik akceptacji będzie bardzo niski, ponieważ propozycje prawdopodobnie trafią do regionów o znacznie mniejszej gęstości prawdopodobieństwa, a więc będą bardzo małe i ponownie łańcuch będzie się bardzo powoli zbiegał. Zazwyczaj dostraja się rozkład propozycji tak, aby algorytmy akceptowały około 30% wszystkich próbek – zgodnie z teoretycznymi szacunkami wspomnianymi w poprzednim akapicie.

Wynik trzech łańcuchów Markowa działających na funkcji 3D Rosenbrocka przy użyciu algorytmu Metropolis-Hastings. Algorytm pobiera próbki z regionów, w których prawdopodobieństwo a posteriori jest wysokie, a łańcuchy zaczynają się mieszać w tych regionach. Przybliżona pozycja maksimum została podświetlona. Czerwone punkty to te, które pozostają po procesie wygrzewania. Wcześniejsze zostały odrzucone.

Zobacz też

Bibliografia

Uwagi

Dalsza lektura