Tekst źródłowy strony Algorytm Metropolisa-Hastingsa

[[Plik:Metropolis algorithm convergence example.png|mały|Coraz dokładniejsze szacowanie rozkładu przy pomocy algorytmu Metropolisa-Hastingsa z użyciem rosnącej liczby powtórzeń]]
[[Plik:MH Rosenbrock movie ND.gif|mały|Znalezienie maksimum funkcji przy pomocy algorytmu Metropolisa-Hastingsa]]
'''Algorytm Metropolisa-Hastingsa''' – metoda statystyczna typu MCMC ([[Próbkowanie Monte Carlo łańcuchami Markowa|próbkowania Monte Carlo łańcuchami Markowa]]), pozwalająca stochastycznie oszacowywać [[Całka|całki]] (takie jak [[estymator]]y) i [[Rozkład prawdopodobieństwa|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 [[Wnioskowanie bayesowskie|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.<ref>{{Cytuj pismo |autor=W.K. Hastings |tytuł=Monte Carlo sampling methods using Markov chains and their applications |czasopismo=Biometrika |wolumin=57 |numer=1 |strony=97–109 |doi=10.1093/biomet/57.1.97}}</ref> [[Metoda Monte Carlo|Metody Monte Carlo]] tego rodzaju stosowali także [[Enrico Fermi]] i [[Stanisław Ulam]] w trakcie tajnych prac nad [[Projekt Manhattan|Projektem Manhattan]]<ref name=":0" />.

Systemy wielowymiarowe obarczone są zjawiskiem ''[[Przekleństwo wymiarowości|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ńcuch Markowa|ł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<ref name=":0">{{Cytuj |autor=Christian P. Robert |tytuł=The Metropolis-Hastings algorithm |czasopismo=arXiv - Computation |doi=10.48550/arXiv.1504.01896 |arxiv=1504.01896 |data=2015-04-08 |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}}</ref>.

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łądzenie losowe|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<ref name=":0" /><ref>{{Cytuj |autor=Michael Betancourt |tytuł=A Conceptual Introduction to Hamiltonian Monte Carlo |arxiv=1701.02434 |doi= 10.48550/arXiv.1701.02434 |czasopismo=arXiv - Methodology |data=2017-01-09}}</ref>.

== Algorytm ==
Algorytm polega na powtarzaniu coraz dokładniejszych iteracji próbkowania rozkładu P(x), w oparciu o błądzenie losowe według wybranej funkcji <math>g,</math> 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 ===
Pierwszy stan <math>x</math> łańcucha Markowa jest wybierany losowo z dziedziny danego rozkładu losowego.

=== 2. Wybranie następnego stanu ''x''′ ===
Kolejny stan jest wybierany według określonej przez użytkownika funkcji błądzenia losowego <math>g(x'|x).</math> Funkcja <math>g</math> może być oparta np. o [[rozkład normalny]] lub [[Rozkład Studenta|rozkład ''t'' Studenta]] (a w przypadku funkcji wielu zmiennych o  [[wielowymiarowy rozkład normalny]] lub wielowymiarowy [[Rozkład Studenta|rozkład <math>t</math> Studenta]]).

=== 3. Przyjęcie lub odrzucenie stanu ''x''′ ===
Ocena nowego stanu odbywa się według wybranego przez użytkownika kryterium, np. Metropolis proponuje funkcję:
<math>A(x' | x) = \min\left(1,\frac{P(x')}{P(x)}\frac{g(x|x')}{g(x'|x)}\right).</math> 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 ===
Testowanie lokalnych stanów <math>x'</math> odbywa się przyjętą przez użytkownika <math>T</math> liczbę powtórzeń.

=== 5. Przyjęcie nowego stanu <math>x</math> i powrót do punktu 2 ===
Po przetestowaniu stanu <math>x'</math> jest on zapisany jako nowy <math>x</math> i algorytm szuka kolejnego punktu <math>x'.</math>

== Przykład próbkowania ==
Poniżej podano kod w języku [[python]] próbkowania [[Rozkład normalny|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 <u>darmowego notatnika</u> [https://colab.google/ colab google] online.
<syntaxhighlight lang="python" line="1">
# 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()
</syntaxhighlight>
== Oprogramowanie dla metod MCMC ==
Popularną implementacją algorytmu Metropolisa-Hastingsa i innych zaawansowanych metod MCMC jest [[Wolne i otwarte oprogramowanie|darmowe i otwarte oprogramowanie]] STAN, dostępne na przykład w otwartym pakiecie statystycznym [[R (język programowania)|R]], i używane do metod wnioskowania bayesowskiego.

== Zobacz też ==

* [[metoda Monte Carlo]]
* [[metoda Eulera-Maruyamy]]
* [[równanie Langevina]]

== Przypisy ==
{{Przypisy}}

{{Kontrola autorytatywna}}

[[Kategoria:Metody numeryczne]]
[[Kategoria:Modelowanie statystyczne]]