[PYTHON] [Statistiques] Je vais expliquer l'échantillonnage par la méthode de Monte Carlo en chaîne de Markov (MCMC) avec animation.
Il s'agit d'un article qui implémente et explique la méthode Monte Carlo par chaîne de Markov en Python.
«Statistiques de calcul II Méthode de Monte Carlo en chaîne de Markov et ses environs» à la p16
La meilleure façon d'avoir une idée du contenu de cette section est dans n'importe quel langage informatique. Essayez de mettre en œuvre ce qui a été dit ici à partir d'une ardoise vierge.
Alors j'ai essayé docilement. Comme c'est un gros problème, je vais vous expliquer le code et son fonctionnement.
Je vais d'abord montrer l'animation et l'intrigue résultantes: kissing_closed_eyes:
(Période de rodage: 1 à 30 [les données de cette période sont tracées avec des couleurs plus claires], jusqu'à 150 échantillons, y compris le rejet)
Tracez le résultat d'un échantillonnage répété 10 000 fois. (Dont, Burn-in: 2000 fois)
introduction
Tout d'abord, importez les bibliothèques requises.
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)
Distribution cible à échantillonner cette fois
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) \\
Utilisez la distribution normale bidimensionnelle moins les constantes de normalisation. Dans cette distribution, les densités de probabilité sont comparées, de sorte que la partie constante peut être ignorée. Définissez une fonction appelée P (・) en Python et utilisez-la comme distribution cible.
#Distribution cible: fonction de densité de probabilité de la distribution normale bidimensionnelle moins constante de normalisation
def P(x1, x2, b):
assert np.abs(b) < 1
return np.exp(-0.5*(x1**2 - 2*b*x1*x2 + x2**2))
Définissez divers paramètres.
# parameters
b = 0.5 #Covariance de la distribution d'objets
delta = 1 #Écart type de la distribution proposée
dist_type = "norm" #Types de distribution proposée("norm" or "unif")
print "Types de distribution proposée:", dist_type
num_frame = 150. #Nombre total d'images pour l'animation
#Liste pour stocker les résultats d'échantillonnage
sample = []
# 1:accept, 0:Liste à stocker rejet
acc_rej = []
Après avoir défini les paramètres, nous commencerons le processus d'échantillonnage et de dessin de l'animation. J'ai mis des commentaires sur les points importants.
#Définition de la position initiale et stockage dans la liste des résultats d'échantillonnage
current = (3, 3)
sample.append(current)
#Une fonction qui dessine chaque image de l'animation
def animate(nframe):
global current, acc_rej
print nframe, #Afficher la progression
#Sélection de l'étape suivante par distribution proposée
# dist_type: "norm":distribution normale/ "unif": Distribution uniforme
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))
#Rapport de densité de probabilité de la distribution cible à chaque position ...[[1]]
P_prev = P(current[0], current[1], b) #Densité de probabilité de la distribution cible à la position actuelle(Valeur numérique proportionnelle à)
P_next = P(next[0], next[1], b) #Densité de probabilité de la distribution cible au prochain poste candidat(Valeur numérique proportionnelle à)
#Prenez le rapport des deux valeurs ci-dessus
r = P_next/P_prev
#Accepter en haut à gauche du graphique/Afficher le cadre pour afficher Rejeter
ax = fig.add_subplot(111)
rect = plt.Rectangle((-3.8,3.2), 1.1, .5,fc="#ffffff", zorder=nframe)
ax.add_patch(rect)
#Tracez une ligne pointillée représentant la trajectoire du mouvement de la position actuelle à la position candidate suivante
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-Lorsque le nombre aléatoire uniforme de 1 est supérieur à r, l'état est mis à jour.
current = copy.copy(next)
#Remplissez la liste avec les valeurs échantillonnées.
sample.append(current)
if nframe < num_frame*.2:
#Les 20 premières itérations%Est brûler-Pensez-y comme une période(La couleur du tracé est grisée)
alpha = 0.2
else:
#Restaure la densité des points pendant la période normale
alpha = 0.8
#Enregistrer accepter
acc_rej.append(1)
#Adopté(Accept)Alors tracez les points.
plt.scatter(current[0], current[1], alpha=alpha)
plt.text(-3.7, 3.35, "Accept", zorder=nframe, fontdict={'color':"b"})
else:
# 0-Si le nombre aléatoire uniforme de 1 est inférieur à r, il est rejeté.
#En cas de rejet, tracez la marque x.
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:
#Enregistrer le rejet
acc_rej.append(0)
if nframe >= num_frame*.2:
plt.title("cnt:{}".format(nframe+1))
else:
plt.title("cnt:{} [burn-in]".format(nframe+1))
#Définition de la plage de dessin du graphique
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)
#Adopté(Accept)Calcul des tarifs
print "Accept ratio:{0:.5f}".format(np.mean(acc_rej))
Le taux d'acceptation était d'environ 60%. Pas très cher: cher:
Cette méthode de Metropolis Hastings sélectionne la prochaine destination de transition en fonction de la distribution proposée à chaque échantillonnage. Ce qui suit montre la distribution proposée pour sélectionner l'échantillon suivant lorsqu'il est exécuté jusqu'au 40e échantillonnage, avec des courbes de niveau en fonction de sa densité de probabilité. Le point rouge est la ** position actuelle **. Le prochain point candidat qui peut être obtenu à partir de la distribution proposée maintenant
x_1^{(t+1)} = x_1^{(t)} + N(0, 1) \\
x_2^{(t+1)} = x_2^{(t)} + N(0, 1)
Par conséquent, de telles lignes de contour peuvent être dessinées. Étant donné que cette distribution proposée génère au hasard une valeur selon $ N (0, 1) $ quelle que soit la position actuelle, par exemple, comme le montre la figure ci-dessous, la prochaine destination de transition est définie dans le sens de la faible densité de probabilité de la distribution cible. Vous pouvez également choisir.
Vous trouverez ci-dessous une image de la section transversale lorsqu'elle est coupée verticalement le long de la ligne rouge de la figure ci-dessus. La ligne bleue est la distribution cible et la ligne verte est la distribution proposée. Maintenant, un nombre aléatoire est généré à partir de la distribution proposée dont le centre est la position actuelle, et le point bleu est sélectionné comme candidat pour la destination de transition.
Maintenant, je voudrais comparer cette «densité de probabilité de la distribution cible à la position actuelle» et «densité de probabilité de la distribution cible à la destination de transition candidate». Les deux lignes roses ci-dessous correspondent à chacune. Intuitivement, les candidats pour cette destination de transition sont des lieux à faible densité de probabilité compte tenu de la distribution cible, on pense donc qu'ils devraient rarement être réalisés. Maintenant que nous utilisons des nombres aléatoires de la distribution proposée, la densité de probabilité de cette distribution cible n'est pas prise en compte. Afin de refléter la densité de probabilité de la distribution cible, nous pensons que ce candidat de destination de transition sera reflété par l'acceptation et le rejet (Rejet) en fonction d'une certaine règle.
Prenons un ratio et comparons. Si le ratio est $ r $
Peut être exprimé comme. Sur le code Python
#Rapport de densité de probabilité de la distribution cible à chaque position ...[[1]]
P_prev = P(current[0], current[1], b) #Densité de probabilité de la distribution cible à la position actuelle(Valeur numérique proportionnelle à)
P_next = P(next[0], next[1], b) #Densité de probabilité de la distribution cible au prochain poste candidat(Valeur numérique proportionnelle à)
#Prenez le rapport des deux valeurs ci-dessus
r = P_next/P_prev
C'est la partie de.
Ce $ r $ prend un nombre supérieur ou égal à 0. La règle d'adoption est déterminée en interprétant ce ratio comme une certaine probabilité d'acceptation.
Premièrement, si $ r \ ge 1 $ est toujours adopté.
Si $ 0 \ le r <1 $, il est traité de sorte que la valeur de $ r $ puisse être considérée comme la probabilité d'acceptation en la comparant au nombre aléatoire uniforme de [0,1].
Dans ce code Python,
if r > 1 or r > rd.uniform(0, 1): #・ ・ ・[[2]]
…
La partie correspond.
En répétant cela, un échantillonnage pour la distribution cible peut être effectué. La preuve théorique que cela fonctionne est la «Méthode de Monte Carlo de la chaîne de Markov de la statistique informatique II et ses environs». Se il vous plaît se référer.
10000 exécutions d'échantillonnage
L'animation ci-dessus était petite avec 150 échantillons, donc échantillonnons 10000 et tracons-la pour que la distribution cible soit plus visible.
# 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-Lorsque le nombre aléatoire uniforme de 1 est supérieur à r, l'état est mis à jour.
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()
La matrice de co-distribution de distribution
[[1, 0.5],
[0.5, 1]]
L'intrigue ressemble à une distribution normale de: rire:
Transition de valeur moyenne
Regardons également la transition de la valeur moyenne en échantillonnant $ x_1 $ et $ x_2 $. Vous pouvez voir que la valeur moyenne s'approche progressivement de 0 comme prévu. Juste au moment où la valeur moyenne devenait 0, elle atteignait 10000 (y compris 2000 burn-ins), il peut donc être correct d'augmenter un peu plus le nombre d'échantillons.
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()
Histogramme des résultats d'échantillonnage
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()
Ensemble de code
Le code utilisé dans cet article peut être trouvé à [ici] sur Github (https://github.com/matsuken92/Qiita_Contents/blob/master/Bayes_chap_04/metropolis-multi-normal.ipynb).
référence
"Statistiques informatiques II Méthode Monte Carlo en chaîne de Markov et ses environs" Partie I "Principes de base de la méthode de Monte Carlo en chaîne de Markov" (Yukito Iba) https://www.iwanami.co.jp/.BOOKS/00/0/0068520.html
Méthode Markov Chain Monte Carlo (MCMC) à comprendre avec visualisation http://d.hatena.ne.jp/hoxo_m/20140911/p1 ⇒ @hoxo_m Oyabun avait déjà écrit une animation similaire sur son blog ... Désolé pour la deuxième décoction ...
Préférences pour générer des GIF animés à partir de Python sur Mac http://qiita.com/kenmatsu4/items/573ca0733b192d919d0e
Recommended Posts