[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) metropolis_norm_1-compressor.gif

Tracez le résultat d'un échantillonnage répété 10 000 fois. (Dont, Burn-in: 2000 fois)

mcmc10000-compressor.png

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))

metropolis_norm_1-compressor.gif

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.

proposal2-compressor.png

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.

prop_target01-compressor.png

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.

prop_target02-compressor.png

Prenons un ratio et comparons. Si le ratio est $ r $

r.png

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()

mcmc10000-compressor.png

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()

mean_2-compressor.png

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()

hist_02-compressor.png

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

[Statistiques] Je vais expliquer l'échantillonnage par la méthode de Monte Carlo en chaîne de Markov (MCMC) avec animation.
La première méthode de Monte Carlo en chaîne de Markov par PyStan
[Statistiques] Visualisez et comprenez la méthode Hamiltonian Monte Carlo avec animation.
Comprendre la méthode Metropolitan Hasting (une des méthodes de la méthode Monte Carlo en chaîne de Markov) avec implémentation
Estimation de π par la méthode de Monte Carlo
Jeu de compression Dominion analysé par la méthode de Monte Carlo
PRML Chapitre 11 Implémentation Python Monte Carlo Chaîne de Markov
Résumé de chacune des méthodes de la chaîne de Markov et de Monte Carlo
Méthode de Monte Carlo
Comparaison de vitesse de chaque langue par la méthode de Monte Carlo
Recherche du rapport de circonférence avec une fonction à 3 lignes [méthode Python / Monte Carlo]