PRML Chapitre 11 Implémentation Python Monte Carlo Chaîne de Markov
Cette fois, nous avons implémenté l'algorithme Metropolis, qui est un exemple typique de la chaîne de Markov Monte Carlo (MCMC). Cette méthode est souvent utilisée lorsque vous souhaitez échantillonner non seulement des distributions de probabilité célèbres telles que la distribution gaussienne et une distribution uniforme, mais également des distributions avec des formes plus complexes.
Algorithme de Metropolis
Considérons l'échantillonnage à partir d'une distribution de probabilité $ p (x) = {1 \ over Z_p} \ tilde {p} (x) $. Vous n'avez pas besoin de connaître la constante de standardisation $ Z_p $. Lors de la recherche de la distribution de probabilité avec le théorème de Bayes, nous ne connaissons souvent pas la constante standardisée.
Vous ne pouvez pas échantillonner directement à partir d'ici car nous ne connaissons que la fonction non standardisée $ \ tilde {p} (\ cdot) $. Par conséquent, nous préparons une autre distribution de probabilité (par exemple, la distribution gaussienne) qui peut être directement échantillonnée, qui est appelée une distribution proposée. Cependant, la distribution proposée est symétrique.
Méthode Flow of Metropolis
- Définissez la valeur initiale $ x_1 $
- Échantillon d'un échantillon de candidat ($ x ^ \ * $) de la distribution proposée centrée sur $ x_i $
- $ \ min (1, {\ tilde {p} (x ^ *) \ over \ tilde {p} (x_i)}) $ avec une probabilité de $ x_ {i + 1} = x ^ \ * $ Si ce n'est pas accepté, définissez $ x_ {i + 1} = x_i $
- Soit la série $ \ {x_n \} _ {n = 1} ^ N $ obtenue en répétant les étapes 2 et 3 un échantillon de la distribution de probabilité $ p (x) $. C'est une méthode très simple.
La série d'échantillons obtenue par cette procédure a une très forte corrélation avec les échantillons avant et après elle. Puisque nous voulons souvent des échantillons indépendants pour l'échantillonnage, ne conserver que tous les M échantillons de la série d'échantillons affaiblit la corrélation.
code
Bibliothèque
Utilisez matplotlib et numpy
import matplotlib.pyplot as plt
import numpy as np
Méthode Metropolis
class Metropolis(object):
def __init__(self, func, ndim, proposal_std=1., sample_rate=1):
#Fonctions non standardisées
self.func = func
#Dimension des données
self.ndim = ndim
#Distribution proposée(Distribution gaussienne)Écart type de
self.proposal_std = proposal_std
#Échantillons dilués
self.sample_rate = sample_rate
#échantillonnage
def __call__(self, sample_size):
#Réglage de la valeur initiale
x = np.zeros(self.ndim)
#Liste pour contenir des échantillons
samples = []
for i in xrange(sample_size * self.sample_rate):
#Échantillonnage à partir de la distribution proposée
x_new = np.random.normal(scale=self.proposal_std, size=self.ndim)
x_new += x
#Formule PRML(11.33)Calculez la probabilité que l'échantillon de candidats soit accepté
accept_prob = self.func(x_new) / self.func(x)
if accept_prob > np.random.uniform():
x = x_new
#Tenir l'échantillon
if i % self.sample_rate == 0:
samples.append(x)
return np.array(samples)
Fonction principale
def main():
#Fonctions non standardisées
def func(x):
return np.exp(-0.5 * np.sum(x ** 2, axis=-1) / 5.)
#Essayez d'abord dans un espace unidimensionnel
print "one dimensional"
#L'écart type de la distribution proposée est de 2 et les échantillons sont éclaircis tous les 10
sampler = Metropolis(func, ndim=1, proposal_std=2., sample_rate=10)
#Échantillonnez 100 pièces en utilisant la méthode Metropolis
samples = sampler(sample_size=100)
#Vérifier la moyenne et la variance de l'échantillon
print "mean", np.mean(samples)
print "var", np.var(samples, ddof=1)
#Exemples de résultats illustrés
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()
#Ensuite, essayez dans un espace bidimensionnel
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()
Code entier
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()
résultat
Les figures ci-dessous montrent les résultats de l'échantillonnage dans l'espace unidimensionnel et bidimensionnel, respectivement.
Les courbes de niveau sont de la fonction de densité de probabilité.
Résultat de sortie du terminal
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]
Dans les deux cas, la moyenne et la variance de l'échantillon sont respectivement proches de la moyenne et de la variance de la population.
À la fin
Il existe d'autres méthodes comme le Hamiltonian Monte Carlo, donc je les implémenterai si j'en ai l'occasion.
Recommended Posts