PRML Kapitel 11 Implementierung der Markov-Kette Monte Carlo Python
Dieses Mal haben wir den Metropolis-Algorithmus implementiert, der ein typisches Beispiel für die Markov-Kette Monte Carlo (MCMC) ist. Diese Methode wird häufig verwendet, wenn Sie nicht nur berühmte Wahrscheinlichkeitsverteilungen wie die Gaußsche Verteilung und die Gleichverteilung, sondern auch Verteilungen mit komplizierteren Formen abtasten möchten.
Metropolis-Algorithmus
Betrachten Sie eine Stichprobe aus einer Wahrscheinlichkeitsverteilung $ p (x) = {1 \ über Z_p} \ tilde {p} (x) $. Sie müssen die Standardisierungskonstante $ Z_p $ nicht kennen. Wenn wir die Wahrscheinlichkeitsverteilung mit dem Bayes-Theorem finden, kennen wir die standardisierte Konstante oft nicht.
Sie können hier nicht direkt probieren, da wir nur die nicht standardisierte Funktion $ \ tilde {p} (\ cdot) $ kennen. Daher bereiten wir eine andere Wahrscheinlichkeitsverteilung (z. B. die Gaußsche Verteilung) vor, die direkt abgetastet werden kann und als vorgeschlagene Verteilung bezeichnet wird. Die vorgeschlagene Verteilung ist jedoch symmetrisch.
Fluss der Metropolenmethode
- Setzen Sie den Anfangswert $ x_1 $
- Probieren Sie einen Beispielkandidaten ($ x ^ \ * $) aus der vorgeschlagenen Verteilung aus, die auf $ x_i $ zentriert ist
- $ \ min (1, {\ tilde {p} (x ^ *) \ over \ tilde {p} (x_i)}) $ mit einer Wahrscheinlichkeit von $ x_ {i + 1} = x ^ \ * $ Wenn nicht akzeptiert, setzen Sie $ x_ {i + 1} = x_i $
- Die durch Wiederholen der Schritte 2 und 3 erhaltene Reihe $ \ {x_n \} _ {n = 1} ^ N $ sei eine Stichprobe aus der Wahrscheinlichkeitsverteilung $ p (x) $. Es ist eine sehr einfache Methode.
Die durch dieses Verfahren erhaltenen Probenserien weisen eine sehr starke Korrelation mit den Proben davor und danach auf. Da wir häufig unabhängige Stichproben für die Stichprobe benötigen, schwächt die Beibehaltung nur aller M Stichproben in der Stichprobenreihe die Korrelation.
Code
Bibliothek
Verwenden Sie Matplotlib und Numpy
import matplotlib.pyplot as plt
import numpy as np
Metropolis-Methode
class Metropolis(object):
def __init__(self, func, ndim, proposal_std=1., sample_rate=1):
#Nicht standardisierte Funktionen
self.func = func
#Datendimension
self.ndim = ndim
#Vorgeschlagene Verteilung(Gaußsche Verteilung)Standardabweichung von
self.proposal_std = proposal_std
#Proben ausdünnen
self.sample_rate = sample_rate
#Probenahme
def __call__(self, sample_size):
#Anfangswerteinstellung
x = np.zeros(self.ndim)
#Liste für Proben
samples = []
for i in xrange(sample_size * self.sample_rate):
#Stichprobe aus der vorgeschlagenen Verteilung
x_new = np.random.normal(scale=self.proposal_std, size=self.ndim)
x_new += x
#PRML-Formel(11.33)Berechnen Sie die Wahrscheinlichkeit, dass Stichprobenkandidaten akzeptiert werden
accept_prob = self.func(x_new) / self.func(x)
if accept_prob > np.random.uniform():
x = x_new
#Probe halten
if i % self.sample_rate == 0:
samples.append(x)
return np.array(samples)
Hauptfunktion
def main():
#Nicht standardisierte Funktionen
def func(x):
return np.exp(-0.5 * np.sum(x ** 2, axis=-1) / 5.)
#Versuchen Sie es zuerst im eindimensionalen Raum
print "one dimensional"
#Die Standardabweichung der vorgeschlagenen Verteilung beträgt 2, und die Proben werden alle 10 verdünnt
sampler = Metropolis(func, ndim=1, proposal_std=2., sample_rate=10)
#Probieren Sie 100 Stück nach der Metropolis-Methode
samples = sampler(sample_size=100)
#Überprüfen Sie den Stichprobenmittelwert und die Varianz
print "mean", np.mean(samples)
print "var", np.var(samples, ddof=1)
#Illustrierte Beispielergebnisse
x = np.linspace(-10, 10, 100)[:, None]
y = func(x) / np.sqrt(2 * np.pi * 5.)
plt.plot(x, y, label="probability density function")
plt.hist(samples, normed=True, alpha=0.5, label="normalized sample histogram")
plt.scatter(samples, np.random.normal(scale=0.001, size=len(samples)), label="samples")
plt.xlim(-10, 10)
plt.show()
#Versuchen Sie es als nächstes im zweidimensionalen Raum
print "\ntwo dimensional"
sampler = Metropolis(func, 2, proposal_std=2., sample_rate=10)
samples = sampler(sample_size=100)
print "mean\n", np.mean(samples, axis=0)
print "covariance\n", np.cov(samples, rowvar=False)
x, y = np.meshgrid(
np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
z = func(np.array([x, y]).reshape(2, -1).T).reshape(100, 100)
plt.contour(x, y, z)
plt.scatter(samples[:, 0], samples[:, 1])
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()
Ganzer Code
mcmc.py
import matplotlib.pyplot as plt
import numpy as np
class Metropolis(object):
def __init__(self, func, ndim, proposal_std=1., sample_rate=1):
self.func = func
self.ndim = ndim
self.proposal_std = proposal_std
self.sample_rate = sample_rate
def __call__(self, sample_size):
x = np.zeros(self.ndim)
samples = []
for i in xrange(sample_size * self.sample_rate):
x_new = np.random.normal(scale=self.proposal_std, size=self.ndim)
x_new += x
accept_prob = self.func(x_new) / self.func(x)
if accept_prob > np.random.uniform():
x = x_new
if i % self.sample_rate == 0:
samples.append(x)
assert len(samples) == sample_size
return np.array(samples)
def main():
def func(x):
return np.exp(-0.5 * np.sum(x ** 2, axis=-1) / 5.)
print "one dimensional"
sampler = Metropolis(func, ndim=1, proposal_std=2., sample_rate=10)
samples = sampler(sample_size=100)
print "mean", np.mean(samples)
print "var", np.var(samples, ddof=1)
x = np.linspace(-10, 10, 100)[:, None]
y = func(x) / np.sqrt(2 * np.pi * 5.)
plt.plot(x, y, label="probability density function")
plt.hist(samples, normed=True, alpha=0.5, label="normalized sample histogram")
plt.scatter(samples, np.random.normal(scale=0.001, size=len(samples)), label="samples")
plt.xlim(-10, 10)
plt.show()
print "\ntwo dimensional"
sampler = Metropolis(func, 2, proposal_std=2., sample_rate=10)
samples = sampler(sample_size=100)
print "mean\n", np.mean(samples, axis=0)
print "covariance\n", np.cov(samples, rowvar=False)
x, y = np.meshgrid(
np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
z = func(np.array([x, y]).reshape(2, -1).T).reshape(100, 100)
plt.contour(x, y, z)
plt.scatter(samples[:, 0], samples[:, 1])
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()
if __name__ == '__main__':
main()
Ergebnis
Die folgenden Abbildungen zeigen die Ergebnisse der Probenahme im eindimensionalen bzw. zweidimensionalen Raum.
Die Konturlinien haben die Wahrscheinlichkeitsdichtefunktion.
Ergebnis der Terminalausgabe
terminal
one dimensional
mean 0.427558835137
var 5.48086205252
two dimensional
mean
[-0.04893427 -0.04494551]
covariance
[[ 5.02950816 -0.02217824]
[-0.02217824 5.43658538]]
[Finished in 1.8s]
In beiden Fällen liegen der Stichprobenmittelwert und die Varianz nahe am Populationsmittelwert bzw. der Populationsvarianz.
Am Ende
Es gibt andere Methoden wie Hamiltonian Monte Carlo, daher werde ich sie implementieren, wenn ich die Gelegenheit dazu habe.
Recommended Posts