[PYTHON] Verstehen Sie die Metropolitan Hasting-Methode (eine der Methoden in der Markov-Ketten-Monte-Carlo-Methode) mit der Implementierung
Einführung
Zurück erläuterte die Markov-Kette bzw. die Monte-Carlo-Methode und fasste die Umrisse der Markov-Ketten-Monte-Carlo-Methode (allgemein als MCMC-Methode bekannt) zusammen. Es gibt mehrere Algorithmen für die MCMC-Methode, aber dieses Mal möchte ich die Metropolitan Hasting-Methode (allgemein als MH-Methode bekannt) zusammenfassen.
Ich habe diesmal gelernt, indem ich mich auf die folgenden Artikel und Bücher bezog.
-
[Bayes-Statistiken aus den Grundlagen: Eine praktische Einführung in die Hamilton-Monte-Carlo-Methode](https://www.amazon.co.jp/%E5%9F%BA%E7%A4%8E%E3%81%8B%E3 % 82% 89% E3% 81% AE% E3% 83% 99% E3% 82% A4% E3% 82% BA% E7% B5% B1% E8% A8% 88% E5% AD% A6-% E3% 83% 8F% E3% 83% 9F% E3% 83% AB% E3% 83% 88% E3% 83% 8B% E3% 82% A2% E3% 83% B3% E3% 83% A2% E3% 83% B3% E3% 83% 86% E3% 82% AB% E3% 83% AB% E3% 83% AD% E6% B3% 95% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E5% AE% 9F% E8% B7% B5% E7% 9A% 84% E5% 85% A5% E9% 96% 80-% E8% B1% 8A% E7% 94% B0-% E7% A7% 80% E6% A8% B9 / dp / 4254122128)
-
Bayes-Statistiken aus dem Basics Round-Reading-Material Kapitel 4 Metropolis Hastings-Methode
Rückblick auf die Markov-Ketten-Monte-Carlo-Methode
Die Markov-Ketten-Monte-Carlo-Methode (allgemein als MCMC-Methode bekannt) ist eine Methode zum Erhalten einer komplexen Verteilung wie der Wahrscheinlichkeitsverteilung in einem mehrdimensionalen Raum aus ** Stichproben basierend auf der Wahrscheinlichkeitsverteilung **.
Bei der Berechnung des erwarteten Werts einer Funktion $ θ_ {eap} $ handelt es sich in der Regel um eine hochdimensionale Funktion für die Datenanalyse. In diesem Fall ist die Berechnung sehr kompliziert und oft schwer zu lösen. Ziehen Sie daher in Betracht, die Zielverteilung zufällig zu ermitteln (Monte-Carlo-Methode). Die Verteilung der Abtastwerte wird durch die Markov-Kette so eingestellt, dass die gewünschte Verteilung erhalten wird. Diese MCMC-Methode verfügt über die folgenden vier Hauptmethoden und -algorithmen.
- ** Metropolis Hasting Methode (diesmal erklärt) **
- Hamiltonsche Monte-Carlo-Methode
- Gibt Probenahme
- Scheibenprobenahme
Dieses Mal werde ich diese hastige Methode der Metropole zusammenfassen.
Was ist die Metropolis Hasting-Methode?
Berücksichtigen Sie bei der MCMC-Methode die unten gezeigten detaillierten Gleichgewichtsbedingungen für beliebige Wahrscheinlichkeitsvariablen $ θ und θ '$.
f(θ'|θ)f(θ) = f(θ|θ')f(θ') \hspace{1cm}Gleichung 1
In diesem Moment,
Also dieser Übergangskern
Die vorgeschlagene Verteilung ist eine Wahrscheinlichkeitsverteilung, die geeignete Kandidaten als Stichprobe aus der Zielverteilung vorschlägt. Diese vorgeschlagene Verteilung erfüllt nicht unbedingt die Gleichheit, da sie angemessen gewählt wird. Zum Beispiel
q(θ|θ')f(θ') > q(θ'|θ)f(θ)\hspace{1cm}Gleichung 2\\
Es sieht hier aus wie Gleichung 2. Die Methode zur Durchführung einer Wahrscheinlichkeitskorrektur zur Erfüllung dieser Formel als detaillierte Gleichgewichtsbedingung wird als ** Metropolis-Hastings-Algorithmusmethode (allgemein als MH-Methode bekannt) ** bezeichnet.
Ursprünglich ist es ein Algorithmus, der auf dem Gebiet der physikalischen Chemie vorgeschlagen wurde. Es wurde auf die Bayes'sche Statistik angewendet. Gleichung 2 zeigt, dass die Wahrscheinlichkeit, sich zu $ θ $ zu bewegen, größer ist als die Wahrscheinlichkeit, sich zu $ θ '$ zu bewegen. Um jedoch mit der Übergangswahrscheinlichkeit zu korrigieren, führen wir Korrekturkoeffizienten $ c $ und $ c '$ mit positiven Vorzeichen ein.
f(θ|θ') = cq(θ|θ')\\
f(θ'|θ) = c'q(θ'|θ)
Daher nach Korrektur
cq(θ|θ')f(θ')=c'q(θ'|θ)f(θ)
Und die Integration wurde hergestellt. Aber auch in diesem Fall
- Es gibt zwei Korrekturkoeffizienten [tex: c, c '] und es ist redundant
- Es ist schwierig, probabilistische Maßnahmen zu ergreifen, da der Korrekturkoeffizient nicht innerhalb von 0 oder mehr und 1 oder weniger liegt.
Es gibt Probleme wie. Deshalb,
r = \frac{c}{c'} = \frac{q(θ'|θ)f(θ)}{q(θ|θ')f(θ')}\\
r' =\frac{c'}{c'} = 1
Es wird sein. Dies ersetzen
rq(θ|θ')f(θ')=r'q(θ'|θ)f(θ)
Weil es wird
MH-Algorithmus
Nun ist der MH-Algorithmus unten.
- Verwenden Sie die vorgeschlagene Verteilung $ q (| θ ^ {t}) $, um eine Zufallszahl $ a $ zu generieren.
- Bestimmen Sie die folgenden Sätze.
q(a|θ^{t})f(θ^{t}) > q(θ^{t}|a)f(a)
- True: Aus der Ungleichungsbedingung von Gleichung 1 wird bestimmt, dass $ θ = a $, $ θ '= θ ^ {t} $. In diesem Fall ist die vorgeschlagene Verteilung $ q (θ | θ ') $, so dass eine Wahrscheinlichkeitskorrektur unter Verwendung des Koeffizienten $ r $ auf der linken Seite von Gleichung 2 durchgeführt werden muss.
r =\frac {q(θ^{t}|a)f(a) }{q(a|θ^{t})f(θ^{t})}
Wird berechnet, akzeptiert $ a $ mit einer Wahrscheinlichkeit von $ r $ und setzt $ θ ^ {t + 1} = a $. Verwerfe $ a $ mit einer Wahrscheinlichkeit von $ 1-r $ und setze $ θ ^ {t + 1} = θ ^ t $.
- False: Aus der Ungleichungsbedingung von Ausdruck 1 wird bestimmt, dass $ θ '= a $, $ θ = θ ^ {t} $. In diesem Fall ist die vorgeschlagene Verteilung $ q (θ '| θ) $, so dass eine Wahrscheinlichkeitskorrektur unter Verwendung des Koeffizienten $ r' $ auf der rechten Seite von Gleichung 2 durchgeführt werden muss. Akzeptiere immer $ a $ mit einer Wahrscheinlichkeit von $ 1 (= r ') $ und setze $ θ ^ {t + 1} = a $.
- Setzen Sie $ t = t + 1 $ und kehren Sie zu 1 zurück.
Zusammenfassend kann angenommen werden, dass ** den vorgeschlagenen Kandidatenpunkt $ a $ mit einer Wahrscheinlichkeit von $ min (1, r) $ $ (θ ^ {t + 1} = a) $ akzeptiert oder an Ort und Stelle $ θ ^ bleibt {t + 1} = θ ^ {t} $ wird beliebig oft wiederholt **.
Der Umriss der Bewegung ist in der folgenden Abbildung dargestellt. $ Θ $ macht es einfach, sich in Richtung der Zielverteilung $ f (θ) $ zu bewegen.
Implementieren und überprüfen (unabhängige MH-Methode)
Dann möchte ich diese MH-Methode tatsächlich implementieren. Es gibt zwei Arten von MH-Methoden: die unabhängige MH-Methode und die Random-Walk-MH-Methode. Zunächst zur unabhängigen MH-Methode. Das Thema ist wie folgt.
Es wird angenommen, dass die hintere Verteilung der Population $ θ $ die Gammaverteilung von $ f ist (θ | α = 11, λ = 13) $. Dann ist die vorgeschlagene Verteilung eine Normalverteilung $ q (θ) $ mit einem Durchschnitt von 1 und einer Varianz von 0,5. Hier ist der Korrekturkoeffizient
r =\frac {q(θ^{t})f(a) }{q(a)f(θ^{t})}
Und. Das Programm ist unten (der Import der Bibliothek usw. ist nicht geschrieben. Bitte beachten Sie den Volltext-Link später).
theta = []
# Initial value
current = 3
theta.append(current)
n_itr = 100000
for i in range(n_itr):
#Zufällige Generierung aus der vorgeschlagenen Verteilung
a = rand_prop()
if a < 0:
theta.append(theta[-1])
continue
r = (q(current)*f_gamma(a)) / (q(a)*f_gamma(current))
assert r > 0
if r < 0:
#reject
theta.append(theta[-1])
continue
if r >= 1 or r > st.uniform.rvs():#Generiere eine Zufallszahl und akzeptiere sie, wenn sie größer als r ist
# Accept
theta.append(a)
current = a
else:
#Reject
theta.append(theta[-1])
Der Punkt ist, wie man ausdrückt, ob man mit einer Wahrscheinlichkeit von $ r $ akzeptiert.
if r >= 1 or r > st.uniform.rvs()
Es wird entschieden, dass die Methode "st.uniform.rvs ()" Zufallszahlen in einer gleichmäßigen Verteilung von $ 0 bis $ 1 generiert und größere akzeptiert.
Überprüfen Sie anhand dieses Programms, ob die posteriore Verteilung reproduziert werden kann.
plt.title("Histgram from MCMC sampling.")
plt.hist(theta[n_burn_in:], bins=30, normed=True)
xx = np.linspace(0, 2.5,501)
plt.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))
plt.show()
Die als Histogramm blau dargestellten Daten sind die Häufigkeitsverteilung von $ θ $, die nach der MH-Methode verschoben wurde. Es wird normalisiert und korrigiert, so dass die Summe 1 ist. Da die orange Farbe die hintere Verteilung ist, die wir ursprünglich reproduzieren wollten, können wir bestätigen, dass sie reproduziert wurde.
Nun, es sieht so aus, als ob es ohne Probleme gemacht wurde, aber in der vorgeschlagenen Verteilung betrug der Durchschnitt 1 und die Varianz 0,5. Dies ist willkürlich und funktioniert möglicherweise nicht.
Wenn beispielsweise eine Normalverteilung mit einem Mittelwert von 3 und einer Varianz von 0,5 die vorgeschlagene Verteilung ist,
Es wird sich so verschieben. Tatsächlich gibt es bei der unabhängigen MH-Methode das Problem, dass die Verteilung nicht gut konvergiert, es sei denn, die Verteilung nahe der Zielverteilung ist eine stetige Verteilung. Mit anderen Worten, die unabhängige MH-Methode macht einen großen Unterschied in den Ergebnissen bis zur Konvergenz in Abhängigkeit von der Qualität der vorgeschlagenen Verteilung. Dieses Mal wurde die posteriore Verteilung der Einfachheit halber geklärt, aber in vielen Fällen ist die Verteilung in der tatsächlichen Datenanalyse unbekannt. Daher müssen wir uns überlegen, wie wir dieses Problem lösen können.
Random Walk MH-Methode
Ein Algorithmus, der als Random-Walk-MH-Methode bezeichnet wird, wurde als Methode zur Lösung dieses Problems vorgeschlagen.
Kandidaten für Vorschläge
a =θ^{t} + e
Wird besorgt. Für $ e $ wird eine Normalverteilung mit einem Durchschnitt von 0 oder eine Gleichverteilung verwendet.
Wenn Sie eine symmetrische Verteilung wie eine Normalverteilung oder eine Gleichverteilung für die vorgeschlagene Verteilung wählen,
q(a|θ^{t}) = q(θ^{t}|a)
Und die vorgeschlagene Verteilung. Daher ist der Korrekturkoeffizient r bei der Random-Walk-MH-Methode
r = \frac {f(a)}{f(θ^{t})}
Und die vorgeschlagene Verteilung verschwindet immer.
Implementieren und überprüfen (Random Walk MH-Methode)
Lassen Sie uns nun das gleiche Problem wie zuvor implementieren und bestätigen.
# MCMC sampling
theta = []
current = 5
theta.append(current)
n_itr = 100000
prop_m = 5
prop_sd = np.sqrt(0.10)
x = np.linspace(-1,5,501)
def rand_prop1(prop_m):
return st.norm.rvs(prop_m, prop_sd)
for i in range(n_itr):
a = rand_prop1(current) #Zufällige Generierung aus der vorgeschlagenen Verteilung
r = f_gamma(a) / f_gamma(current)
if a < 0 or r < 0:
#reject
theta.append(theta[-1])
pass
elif r >= 1 or r > st.uniform.rvs():
# Accept
theta.append(a)
current = a
status = "acc"
col = "r"
else:
#Reject
theta.append(theta[-1])
pass
Die Grafikbeschreibung ist unten.
plt.title("Histgram from MCMC sampling.")
n, b, p = plt.hist(theta[n_burn_in:], bins=50, density=True)
xx = np.linspace(0, 2.5,501)
plt.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))
plt.show()
kk = 11.
tt = 13.
print ("sample mean:{0:.5f}, sample std:{1:.5f}".format(np.mean(theta), np.std(theta)))
#Theoretischer Erwartungswert / Dispersion
print("mean:{0:.5f}, std:{1:.5f}".format( kk/tt, np.sqrt(kk*((1/tt)**2))) )
Es stellte sich heraus, dass die Verteilung ordentlich reproduziert wurde. Im Vergleich zu der zuvor erwähnten unabhängigen MH-Methode stellte ich fest, dass die Verwendung dieses Algorithmus für die MH-Methode mit zufälliger Wanderung sicher zu sein scheint.
Das Problem der MH-Methode ist jedoch
Ich möchte sagen, dass dies eine Erleichterung ist, aber in Wirklichkeit sind noch Probleme zu lösen. Schließlich müssen Sie den Wert von $ e $ richtig angeben, daher ist auch eine hyperparameterische Denkweise erforderlich. Wenn die Varianz klein ist, ist der Prozentsatz der akzeptierten Übergänge hoch, aber das Fortschreiten im Zustandsraum ist ein langsamer Zufallslauf, der eine lange Konvergenzzeit erfordert. Andererseits ist die Rückweisungsrate umso höher, je höher die Dispersion ist. Die Hamiltonsche Monte-Carlo-Methode ist ein Weg, dies zu lösen. Ich hoffe, dass das nächste Mal zusammengefasst wird.
Am Ende
Dieses Mal haben wir die Metropolis Hasting-Methode zusammengefasst. Das Überlegen, wie man gegen den Mittelwert der posterioren Verteilung konvergiert, ist wie das Optimieren von Gewichtsparametern in einem neuronalen Netzwerk.
Obwohl die Algorithmen unterschiedlich sind, ist der Zweck, den ich tun möchte, ähnlich. Ich dachte, wenn ich das verstehen könnte, könnte ich den Spaß und die Breite von Statistik und maschinellem Lernen sehen.
Der vollständige Text des Programms wird hier gespeichert. https://github.com/Fumio-eisan/Montecarlo20200516