[PYTHON] [Statistik] Ich werde die Abtastung nach der Markov-Ketten-Monte-Carlo-Methode (MCMC) mit Animation erläutern.
Dies ist ein Artikel, der die Markov-Ketten-Monte-Carlo-Methode in Python implementiert und erklärt.
"Berechnungsstatistik II Markov-Ketten-Monte-Carlo-Methode und ihre Umgebung" auf S. 16
Der beste Weg, um ein Gefühl für den Inhalt dieses Abschnitts zu bekommen, ist in jeder Computersprache. Versuchen Sie, das hier Gesagte von einer leeren Tafel aus umzusetzen.
Also habe ich es gehorsam versucht. Da es eine große Sache ist, erkläre ich den Code und wie er funktioniert.
Ich werde zuerst die resultierende Animation und Handlung zeigen: kissing_closed_eyes:
(Einbrennzeitraum: 1-30 [Daten in diesem Zeitraum sind mit helleren Farben dargestellt], bis zu 150 Proben einschließlich Ablehnung)
Zeichnen Sie das Ergebnis der 10.000-mal wiederholten Probenahme auf. (Davon Einbrennen: 2.000 Mal)
Einführung
Importieren Sie zunächst die erforderlichen Bibliotheken.
import numpy as np
import numpy.random as rd
import pandas as pd
import scipy.stats as st
import copy
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import animation as ani
import matplotlib.cm as cm
import seaborn as sns
sns.set(style="whitegrid", palette="muted", color_codes=True)
Zielverteilung, die dieses Mal abgetastet werden soll
f(x,y) = {\sqrt{1-b^2} \over 2\pi } \exp \left( -{1 \over 2} (x_1^2 - 2bx_1x_2 + x_2^2) \right) \\
\propto \exp \left( -{1 \over 2} (x_1^2 - 2bx_1x_2 + x_2^2) \right) \\
Verwenden Sie die zweidimensionale Normalverteilung abzüglich der Normalisierungskonstanten. Bei dieser Verteilung werden die Wahrscheinlichkeitsdichten verglichen, sodass der konstante Teil ignoriert werden kann. Definieren Sie in Python eine Funktion namens P (・) und verwenden Sie diese als Zielverteilung.
#Zielverteilung: Wahrscheinlichkeitsdichtefunktion der zweidimensionalen Normalverteilung minus Normalisierungskonstante
def P(x1, x2, b):
assert np.abs(b) < 1
return np.exp(-0.5*(x1**2 - 2*b*x1*x2 + x2**2))
Definieren Sie verschiedene Parameter.
# parameters
b = 0.5 #Kovarianz der Objektverteilung
delta = 1 #Standardabweichung der vorgeschlagenen Verteilung
dist_type = "norm" #Arten der vorgeschlagenen Verteilung("norm" or "unif")
print "Arten der vorgeschlagenen Verteilung:", dist_type
num_frame = 150. #Gesamtzahl der Bilder für die Animation
#Liste zum Speichern der Stichprobenergebnisse
sample = []
# 1:accept, 0:Liste zum Speichern von Ablehnungen
acc_rej = []
Nachdem Sie die Parameter eingestellt haben, beginnen wir mit dem Abtasten und Zeichnen der Animation. Ich habe Kommentare in die wichtigen Punkte eingefügt.
#Einstellen der Ausgangsposition und Speichern in der Ergebnisliste der Probenahme
current = (3, 3)
sample.append(current)
#Eine Funktion, die jeden Frame der Animation zeichnet
def animate(nframe):
global current, acc_rej
print nframe, #Fortschritt anzeigen
#Auswahl des nächsten Schritts durch vorgeschlagene Verteilung
# dist_type: "norm":Normalverteilung/ "unif": Gleichmäßige Verteilung
if dist_type == "norm":
next = (current[0] + rd.normal(0, delta), current[1] + rd.normal(0, delta))
else:
next = (current[0] + rd.uniform(-delta, delta), current[1] + rd.uniform(-delta, delta))
#Verhältnis der Wahrscheinlichkeitsdichte der Zielverteilung an jeder Position ...[[1]]
P_prev = P(current[0], current[1], b) #Wahrscheinlichkeitsdichte der Zielverteilung an der aktuellen Position(Zahlenwert proportional zu)
P_next = P(next[0], next[1], b) #Wahrscheinlichkeitsdichte der Zielverteilung an der nächsten Kandidatenposition(Zahlenwert proportional zu)
#Nehmen Sie das Verhältnis der beiden oben genannten Werte
r = P_next/P_prev
#Akzeptieren Sie oben links im Diagramm/Zeigen Sie den Rahmen an, um Ablehnen anzuzeigen
ax = fig.add_subplot(111)
rect = plt.Rectangle((-3.8,3.2), 1.1, .5,fc="#ffffff", zorder=nframe)
ax.add_patch(rect)
#Zeichnen Sie eine gepunktete Linie, die den Bewegungspfad von der aktuellen Position zur nächsten Kandidatenposition darstellt
plt.plot([current[0], next[0]], [current[1], next[1]], "k--", lw=.3, color="gray")
if r > 1 or r > rd.uniform(0, 1): #・ ・ ・[[2]]
# 0-Wenn die einheitliche Zufallszahl 1 größer als r ist, wird der Zustand aktualisiert.
current = copy.copy(next)
#Füllen Sie die Liste mit den abgetasteten Werten.
sample.append(current)
if nframe < num_frame*.2:
#Die ersten 20 Iterationen%Ist brennen-Betrachten Sie es als eine In-Periode(Die Farbe des Diagramms ist abgeblendet.
alpha = 0.2
else:
#Stellt die Punktdichte während des normalen Zeitraums wieder her
alpha = 0.8
#Aufnahme akzeptieren
acc_rej.append(1)
#Angenommen(Accept)Zeichnen Sie also die Punkte.
plt.scatter(current[0], current[1], alpha=alpha)
plt.text(-3.7, 3.35, "Accept", zorder=nframe, fontdict={'color':"b"})
else:
# 0-Wenn die einheitliche Zufallszahl 1 kleiner als r ist, wird sie zurückgewiesen.
#Wenn abgelehnt, zeichnen Sie die x-Markierung.
plt.scatter(next[0], next[1], alpha=0.5, color="r", marker="x")
plt.text(-3.7, 3.35, "Reject", zorder=nframe, fontdict={'color':"r"})
if nframe <= num_frame*.2:
#Rekord ablehnen
acc_rej.append(0)
if nframe >= num_frame*.2:
plt.title("cnt:{}".format(nframe+1))
else:
plt.title("cnt:{} [burn-in]".format(nframe+1))
#Festlegen des Zeichenbereichs des Diagramms
plt.xlim(-4, 4)
plt.ylim(-4, 4)
fig = plt.figure(figsize=(8,7))
anim = ani.FuncAnimation(fig, animate, frames=int(num_frame), blit=True)
anim.save('metropolis_norm.gif', writer='imagemagick', fps=3, dpi=96)
#Angenommen(Accept)Ratenberechnung
print "Accept ratio:{0:.5f}".format(np.mean(acc_rej))
Die Akzeptanzquote lag bei ca. 60%. Nicht sehr teuer: teuer:
Diese Metropolis Hastings-Methode wählt bei jeder Abtastung das nächste Übergangsziel gemäß der vorgeschlagenen Verteilung aus. Das Folgende zeigt die vorgeschlagene Verteilung zum Auswählen der nächsten Stichprobe, wenn diese bis zur 40. Abtastung ausgeführt wird, mit Konturlinien entsprechend ihrer Wahrscheinlichkeitsdichte. Der rote Punkt ist die ** aktuelle Position **. Der nächste Kandidatenpunkt, der jetzt aus der vorgeschlagenen Verteilung erhalten werden kann
x_1^{(t+1)} = x_1^{(t)} + N(0, 1) \\
x_2^{(t+1)} = x_2^{(t)} + N(0, 1)
Daher können solche Konturlinien gezeichnet werden. Da diese vorgeschlagene Verteilung zufällig einen Wert gemäß $ N (0, 1) $ unabhängig von der aktuellen Position erzeugt, wie beispielsweise in der folgenden Abbildung gezeigt, wird das nächste Übergangsziel in Richtung der Dichte mit geringer Wahrscheinlichkeit der Zielverteilung festgelegt. Sie können auch wählen.
Unten sehen Sie ein Bild des Querschnitts, wenn er in der obigen Abbildung vertikal entlang der roten Linie geschnitten wird. Die blaue Linie ist die Zielverteilung und die grüne Linie ist die vorgeschlagene Verteilung. Nun wird aus der vorgeschlagenen Verteilung eine Zufallszahl generiert, deren Zentrum die aktuelle Position ist, und der blaue Punkt wird als Kandidat für das Übergangsziel ausgewählt.
Nun möchte ich diese "Wahrscheinlichkeitsdichte der Zielverteilung an der aktuellen Position" und "Wahrscheinlichkeitsdichte der Zielverteilung am Übergangszielkandidaten" vergleichen. Die beiden rosa Linien unten entsprechen jeweils. Intuitiv sind die Kandidaten für dieses Übergangsziel Orte mit geringer Wahrscheinlichkeitsdichte im Lichte der Zielverteilung, so dass angenommen wird, dass sie selten realisiert werden sollten. Da wir nun Zufallszahlen aus der vorgeschlagenen Verteilung verwenden, wird die Wahrscheinlichkeitsdichte dieser Zielverteilung nicht berücksichtigt. Um die Wahrscheinlichkeitsdichte der Zielverteilung widerzuspiegeln, denken wir, dass dieser Übergangszielkandidat durch Akzeptieren und Ablehnen (Rejet) basierend auf einer bestimmten Regel reflektiert wird.
Nehmen wir ein Verhältnis und vergleichen. Wenn das Verhältnis $ r $ ist
Kann ausgedrückt werden als. Auf dem Python-Code
#Verhältnis der Wahrscheinlichkeitsdichte der Zielverteilung an jeder Position ...[[1]]
P_prev = P(current[0], current[1], b) #Wahrscheinlichkeitsdichte der Zielverteilung an der aktuellen Position(Zahlenwert proportional zu)
P_next = P(next[0], next[1], b) #Wahrscheinlichkeitsdichte der Zielverteilung an der nächsten Kandidatenposition(Zahlenwert proportional zu)
#Nehmen Sie das Verhältnis der beiden oben genannten Werte
r = P_next/P_prev
Es ist der Teil von.
Dieses $ r $ nimmt eine Zahl größer oder gleich 0 an. Die Adoptionsregel wird bestimmt, indem dieses Verhältnis als eine bestimmte Akzeptanzwahrscheinlichkeit interpretiert wird.
Erstens, wenn $ r \ ge 1 $ immer übernommen wird.
Wenn $ 0 \ le r <1 $ ist, wird der Wert von $ r $ so verarbeitet, dass er durch Vergleich mit der einheitlichen Zufallszahl von [0,1] als Akzeptanzwahrscheinlichkeit angesehen werden kann.
In diesem Python-Code
if r > 1 or r > rd.uniform(0, 1): #・ ・ ・[[2]]
…
Der Teil entspricht.
Durch Wiederholen dieses Vorgangs kann eine Abtastung für die Zielverteilung durchgeführt werden. Der theoretische Beweis dafür ist "Computational Statistics II Markov-Ketten-Monte-Carlo-Methode und ihre Umgebung". Siehe.
10.000 Stichprobenausführungen
Die obige Animation war mit 150 Samples klein. Probieren Sie also 10000 aus und zeichnen Sie sie so, dass die Zielverteilung besser sichtbar ist.
# parameters
b = 0.5
delta = 1
dist_type = "norm" # "norm" # "unif"
n_samples = 10000
# result
sample = []
#initial state
current = (5, 5)
sample.append(current)
print dist_type
cnt = 1
while cnt < n_samples:
# candidate of next step
if dist_type == "norm":
next = (current[0] + rd.normal(0, delta), current[1] + rd.normal(0, delta))
else:
next = (current[0] + rd.uniform(-delta, delta), current[1] + rd.uniform(-delta, delta))
P_prev = P(current[0], current[1], b)
P_next = P(next[0], next[1], b)
r = P_next/P_prev
if r > 1 or r > rd.uniform(0, 1):
# 0-Wenn die einheitliche Zufallszahl 1 größer als r ist, wird der Zustand aktualisiert.
current = copy.copy(next)
sample.append(current)
cnt += 1
sample = np.array(sample)
plt.figure(figsize=(9,6))
plt.scatter(sample[int(len(sample)*0.2):,0], sample[int(len(sample)*0.2):,1], alpha=0.2)
plt.title("Scatter plot of 2-dim normal random variable with MCMC.")
plt.show()
Die Verteilungs-Co-Verteilungsmatrix
[[1, 0.5],
[0.5, 1]]
Die Handlung sieht aus wie eine Normalverteilung von: Lachen:
Übergang des Durchschnittswertes
Schauen wir uns auch den Übergang des Durchschnittswerts an, während wir $ x_1 $ und $ x_2 $ abtasten. Sie können sehen, dass sich der Durchschnittswert wie erwartet allmählich 0 nähert. Gerade als der Durchschnittswert 0 wurde, erreichte er 10000 (einschließlich 2000 Einbrenner), so dass es in Ordnung sein kann, die Anzahl der Proben etwas weiter zu erhöhen.
ave = [[],[]]
start = len(sample) * 0.2
for i, d in enumerate(np.array(sample[int(start):])):
#print d
for j in range(2):
if i == 0:
ave[j].append(float(d[j]))
else:
ave[j].append( (ave[j][i-1]*i + d[j])/float(i+1) )
plt.figure(figsize=(15, 5))
plt.xlim(0, len(sample[int(start):]))
plt.plot(np.array(ave).T, lw=1)
plt.title("Sequence of x's and y's mean.")
plt.show()
Histogramm des Stichprobenergebnisses
fig = plt.figure(figsize=(15,6))
ax = fig.add_subplot(121)
plt.hist(sample[start:,0], bins=30)
plt.title("x axis")
ax = fig.add_subplot(122)
plt.hist(sample[start:,1], bins=30, color="g")
plt.title("y axis")
plt.show()
Satz von Code
Der in diesem Artikel verwendete Code finden Sie unter [hier] auf Github (https://github.com/matsuken92/Qiita_Contents/blob/master/Bayes_chap_04/metropolis-multi-normal.ipynb).
Referenz
"Computerstatistik II Markov-Ketten-Monte-Carlo-Methode und Umgebung" Teil I "Grundlagen der Markov-Ketten-Monte-Carlo-Methode" (Yukito Iba) https://www.iwanami.co.jp/.BOOKS/00/0/0068520.html
Markov-Ketten-Monte-Carlo-Methode (MCMC) zum Verständnis der Visualisierung http://d.hatena.ne.jp/hoxo_m/20140911/p1 ⇒ @hoxo_m Oyabun hatte bereits eine ähnliche Animation in seinem Blog geschrieben ... Entschuldigung für den zweiten Sud ...
Einstellungen zum Generieren animierter GIFs aus Python auf dem Mac http://qiita.com/kenmatsu4/items/573ca0733b192d919d0e
Recommended Posts