Algorytm Metropolisa-Hastingsa

Algorytm Metropolisa-Hastingsa – metoda statystyczna typu MCMC (próbkowania Monte Carlo łańcuchami Markowa), pozwalająca stochastycznie oszacowywać całki (takie jak estymatory) i rozkłady prawdopodobieństwa dla złożonych systemów, które są zbyt trudne do modelowania analitycznego, np. układów wielowymiarowych. Narzędzie to jest wykorzystywane między innymi we wnioskowaniu bayesowskim i modelowaniu systemów fizycznych. Procedurę opisał po raz pierwszy publicznie zespół Metropolisa w 1953 r., a rozwinął ją Hastings w 1970 r.^[1] Metody Monte Carlo tego rodzaju stosowali także Enrico Fermi i Stanisław Ulam w trakcie tajnych prac nad Projektem Manhattan^[2].

Systemy wielowymiarowe obarczone są zjawiskiem przekleństwa wymiarowości, polegającego na tym, że wraz ze wzrostem liczby wymiarów problemu, liczba obserwacji potrzebnych do oszacowania ich cech rośnie wykładniczo. Metody Monte Carlo wykorzystujące łańcuchy Markowa są w dużej mierze odporne na ten problem, ponieważ nie wymagają rozwiązania analitycznego, ani nie przeszukują całej przestrzeni problemu, lecz posługują się iterowanym próbkowaniem stochastycznym, które z każdą kolejną iteracją coraz dokładniej skupiają się na centralnych obszarach rozkładów^[2].

Wprowadzenie metod MCMC łącznie z rosnącą dostępnością komputerów umożliwiło przełamanie wielu ograniczeń stojących przed analizą złożonych problemów w naukach empirycznych, opartych wcześniej głównie o uproszczone metody ortodoksyjnej statystyki. Dokładność dostępnych oszacowań zależy w nowych metodach od dopasowania parametrów łańcucha Markowa do oczekiwanego rozkładu, oraz od zastosowanej liczby powtórzeń próbkowania.

Algorytm Metropolisa-Hastingsa opisuje jeden z prostych w implementacji sposobów skonstruowania łańcucha Markowa przy pomocy błądzenia losowego, pozwalający na rozwiązanie wielu typowych problemów. W późniejszych latach opisano także metody dopasowane do szerszych klas problemów, np. śledzące kształt rozkładu z wykorzystaniem wektorów pędu, w szczególności hamiltonianów^[2]^[3].

Algorytm[edytuj | edytuj kod]

Algorytm polega 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 podanymi niżej krokami. Algorytm wymaga, aby próbkowany rozkład był stabilny i nieokresowy, w przeciwnym razie łańcuch Markowa może ulec zapętleniu.

1. Losowa inicjalizacja[edytuj | edytuj kod]

Pierwszy stan $x$ łańcucha Markowa jest wybierany losowo z dziedziny danego rozkładu losowego.

2. Wybranie następnego stanu x′[edytuj | edytuj kod]

Kolejny stan jest wybierany według określonej przez użytkownika funkcji błądzenia losowego $g(x'|x).$ Funkcja $g$ może być oparta np. o rozkład normalny lub rozkład t Studenta (a w przypadku funkcji wielu zmiennych o wielowymiarowy rozkład normalny lub wielowymiarowy rozkład $t$ Studenta).

3. Przyjęcie lub odrzucenie stanu x′[edytuj | edytuj kod]

Ocena nowego stanu odbywa się według wybranego przez użytkownika kryterium, np. Metropolis proponuje funkcję: $A(x'|x)=\min \left(1,{\frac {P(x')}{P(x)}}{\frac {g(x|x')}{g(x'|x)}}\right).$ Akceptuje się przejście do nowego położenia x', jeśli proponowane przez nią parametry pozwalają na oszacowanie rozkładu z wyższym prawdopodobieństwem.

4. Powrót do punktu 2[edytuj | edytuj kod]

Testowanie lokalnych stanów $x'$ odbywa się przyjętą przez użytkownika $T$ liczbę powtórzeń.

5. Przyjęcie nowego stanu $x$ i powrót do punktu 2[edytuj | edytuj kod]

Po przetestowaniu stanu $x'$ jest on zapisany jako nowy $x$ i algorytm szuka kolejnego punktu $x'.$

Przykład próbkowania[edytuj | edytuj kod]

Poniżej podano kod w języku python próbkowania rozkładu normalnego, nieunormowanego do 1, za pomocą metody Metropolisa. Dodatkowo dodano wykresy, porównując histogram otrzymany z próbkowania rozkładu z unormowanym rozkładem normalnym. Program można uruchomić, korzystając np. z darmowego notatnika colab google online.

# Metoda Metropolisa
# Rozkład próbkowany - to nieunormowany rozkład normalny. 
# Celem całkowania tego rozkładu jest wyznaczenie stałej normującej go do 1

import numpy as np
import matplotlib.pyplot as plt

mean = 0     # Średnia próbkowanego rozkładu
std_dev = 1  # Odchylenie standardowe próbkowanego rozkładu

def rozkład(x):
    return np.exp(-0.5*((x - mean)/std_dev)**2)

# Rozkład docelowy - to rozkład normalny, unormowany do 1;
# jego wykres wykorzystamy do porównania go z rozkładem próbek (histogram),
# otrzymanych za pomocą metody Metropolisa
def rozkład_norm(x):
    return np.exp(-0.5*((x - mean)/std_dev)**2)/(std_dev*np.sqrt(2*np.pi))

# Liczba próbek
n = 1000000

# Punkt startowy - dowolnie wybramy z dziedziny próbkowanej funkcji rozkład(x)
x_aktualny = -2   

# Inicjalizacja licznika akceptowanych próbek
zakceptowane_probki = 0

# Lista do przechowywania próbek
probki = []

# Główna pętla algorytmu Metropolisa-Hastingsa
for _ in range(n):
    # Generowanie x_proponowany z rozkładu normalnego
    x_proponowany = x_aktualny + np.random.normal(0, 0.02)  
    wspolczynnik_akceptacji = rozkład(x_proponowany) / rozkład(x_aktualny)

    # Losowa decyzja o akceptacji lub odrzuceniu propozycji
    if np.random.uniform(0, 1) <= wspolczynnik_akceptacji:
         x_aktualny = x_proponowany
         zakceptowane_probki += 1
         probki.append(x_aktualny)

# Obliczenie przybliżone całki:
Calka =  zakceptowane_probki / n
print("Przybliżona wartość całki:", Calka)

# ---------------  DODATEK: WYKRESY ---------------
# Dzielenie listy próbek przez Całka, w celu unormowania:
# a) tworzę pustą listę na próbki unormowane przez całkę b) dzielę każdy element
# listy 'probki' przez 'Calka' c) dodaję do listy 'probki_unormowane'
probki_normowane = []
for element in probki:
    probki_normowane.append(element / Calka)

# Wykresy: histogram unormowany próbek i wykres rozkładu normalnego
plt.figure(figsize=(8, 6))
plt.hist(probki_normowane, bins=150, density=True, alpha=0.6, color='blue', label='Próbki')
x_zakres = np.linspace(-5, 5, 1000)
plt.plot(x_zakres, rozkład_norm(x_zakres), color='red', label='Rozkład docelowy')
plt.title('Unormowany histogram próbek z nieunormowanego rozkładu normalnego')
plt.xlabel('x')
plt.ylabel('Prawdopodobieństwo')
plt.legend()
plt.grid(True)
plt.show()

Oprogramowanie dla metod MCMC[edytuj | edytuj kod]

Popularną implementacją algorytmu Metropolisa-Hastingsa i innych zaawansowanych metod MCMC jest darmowe i otwarte oprogramowanie STAN, dostępne na przykład w otwartym pakiecie statystycznym R, i używane do metod wnioskowania bayesowskiego.

Zobacz też[edytuj | edytuj kod]

Przypisy[edytuj | edytuj kod]

↑ W.K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. „Biometrika”. 57 (1). s. 97–109. DOI: 10.1093/biomet/57.1.97.
↑ ^a ^b ^c Christian P.Ch.P. Robert Christian P.Ch.P., The Metropolis-Hastings algorithm, „arXiv - Computation”, 2015, DOI: 10.48550/arXiv.1504.01896, arXiv:1504.01896, Cytat: This Metropolis algorithm, while used in physics, was only generalized by Hastings (1970) and Peskun (1973, 1981) towards statistical applications, as a method apt to overcome the curse of dimensionality penalising regular Monte Carlo methods .
↑ MichaelM. Betancourt MichaelM., A Conceptual Introduction to Hamiltonian Monte Carlo, „arXiv - Methodology”, 2017, DOI: 10.48550/arXiv.1701.02434, arXiv:1701.02434 .

[1] W.K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. „Biometrika”. 57 (1). s. 97–109. DOI: 10.1093/biomet/57.1.97.

[:0-2] Christian P.Ch.P. Robert Christian P.Ch.P., The Metropolis-Hastings algorithm, „arXiv - Computation”, 2015, DOI: 10.48550/arXiv.1504.01896, arXiv:1504.01896, Cytat: This Metropolis algorithm, while used in physics, was only generalized by Hastings (1970) and Peskun (1973, 1981) towards statistical applications, as a method apt to overcome the curse of dimensionality penalising regular Monte Carlo methods .

[3] MichaelM. Betancourt MichaelM., A Conceptual Introduction to Hamiltonian Monte Carlo, „arXiv - Methodology”, 2017, DOI: 10.48550/arXiv.1701.02434, arXiv:1701.02434 .

[1]

[2]

[3]