Comprendre la méthode Metropolitan Hasting (une des méthodes de la méthode Monte Carlo en chaîne de Markov) avec implémentation

introduction

Précédent a expliqué la chaîne de Markov et la méthode de Monte Carlo, respectivement, et a résumé les grandes lignes de la méthode de Monte Carlo par chaîne de Markov (communément appelée méthode MCMC). Il existe plusieurs algorithmes pour la méthode MCMC, mais cette fois, je voudrais résumer la méthode métropolitaine de précipitation (communément appelée méthode MH).

J'ai appris cette fois en me référant aux articles et livres suivants.

Retour sur la méthode Monte Carlo en chaîne de Markov

La méthode de Monte Carlo en chaîne de Markov (communément appelée méthode MCMC) est une méthode permettant d'obtenir une distribution complexe telle que la distribution de probabilité dans un espace multidimensionnel à partir d'un ** échantillonnage basé sur la distribution de probabilité **.  image.png

Lors du calcul de la valeur attendue d'une fonction $ θ_ {eap} $, elle a tendance à être une fonction de grande dimension pour l'analyse des données. Dans ce cas, le calcul est très compliqué et souvent difficile à résoudre, envisagez donc un échantillonnage aléatoire de la distribution cible (méthode de Monte Carlo). La distribution des valeurs échantillonnées est ajustée par la chaîne de Markov pour obtenir la distribution souhaitée.   Cette méthode MCMC comprend les quatre principales méthodes et algorithmes suivants.

Cette fois, je résumerai cette méthode de précipitation de métropole.

Qu'est-ce que la méthode de hachage Metropolis?

Dans la méthode MCMC, considérons les conditions d'équilibre détaillées présentées ci-dessous pour les variables de probabilité arbitraires $ θ et θ '$.

f(θ'|θ)f(θ) = f(θ|θ')f(θ') \hspace{1cm}Équation 1

Cette fois,f(θ'|θ)aussi bien quef(θ|θ')MaisNoyau de transitionC'est une distribution appelée. Distribution postérieure généralement connuef(θ)に対して式1を満たすようなNoyau de transitionをいきなり見つけることは非常に困難です。

Donc ce noyau de transitionf(|)Au lieu deAnalyste appropriéNoyau de transitionq(|)ÀDistribution proposéeJe vais l'utiliser comme.この適当に選ぶことが実は難しさÀ表していることが後の実装で分かります。

La distribution proposée est une distribution de probabilité qui propose des candidats appropriés comme échantillon de la distribution cible. Cette distribution proposée ne satisfait pas nécessairement l'égalité car elle est choisie de manière appropriée. Par exemple

q(θ|θ')f(θ') > q(θ'|θ)f(θ)\hspace{1cm}Équation 2\\

Cela ressemble à l'équation 2 ici. La méthode de correction de probabilité pour satisfaire cette formule en tant que condition d'équilibre détaillée est appelée la ** méthode de l'algorithme Metropolis-Hastings (communément appelée méthode MH) **.

A l'origine, c'est un algorithme proposé dans le domaine de la chimie physique. Il a été appliqué aux statistiques bayésiennes.   L'équation 2 montre que la probabilité de passer à $ θ $ est supérieure à la probabilité de passer à $ θ '$. Cependant, afin de corriger avec la probabilité de transition, nous introduisons des coefficients de correction $ c $ et $ c '$ avec des signes positifs.

f(θ|θ') = cq(θ|θ')\\
f(θ'|θ) = c'q(θ'|θ)

Par conséquent, après correction,

 cq(θ|θ')f(θ')=c'q(θ'|θ)f(θ)

Et l'intégration était établie. Mais même dans ce cas

Il y a des problèmes tels que. Donc,

r = \frac{c}{c'} = \frac{q(θ'|θ)f(θ)}{q(θ|θ')f(θ')}\\
r' =\frac{c'}{c'} = 1

Ce sera. Remplacer ceci

 rq(θ|θ')f(θ')=r'q(θ'|θ)f(θ)

Parce que ça devientrq(θ|θ')Quandr'q(θ'|θ)を遷移核Quandすれば、詳細釣り合い条件を満たすこQuandQuandなります。

Algorithme MH

Maintenant, l'algorithme MH est ci-dessous.

  1. Utilisez la distribution proposée $ q (| θ ^ {t}) $ pour générer un nombre aléatoire $ a $.
  2. Déterminez les propositions suivantes.
q(a|θ^{t})f(θ^{t}) > q(θ^{t}|a)f(a)
r  =\frac {q(θ^{t}|a)f(a) }{q(a|θ^{t})f(θ^{t})}

Est calculé, accepte $ a $ avec une probabilité de $ r $ et définit $ θ ^ {t + 1} = a $. Supprimez $ a $ avec une probabilité de $ 1-r $ et définissez $ θ ^ {t + 1} = θ ^ t $.

  1. Réglez $ t = t + 1 $ et revenez à 1.

En résumé, ** acceptez le point candidat proposé $ a $ avec une probabilité de $ min (1, r) $ $ (θ ^ {t + 1} = a) $, ou restez en place $ θ ^ {t + 1} = θ ^ {t} $ est répété n'importe quel nombre de fois **. Le contour du mouvement est illustré dans la figure ci-dessous. $ Θ $ facilite le passage à la distribution cible $ f (θ) $.   image.png

Mettre en œuvre et vérifier (méthode MH indépendante)

Ensuite, je voudrais mettre en œuvre cette méthode MH. Il existe deux types de méthodes MH: la méthode MH indépendante et la méthode MH à marche aléatoire. Tout d'abord, à propos de la méthode indépendante MH. Le thème est le suivant.

On suppose que la distribution a posteriori de la population $ θ $ est la distribution gamma de $ f (θ | α = 11, λ = 13) $. Ensuite, la distribution proposée est une distribution normale $ q (θ) $ avec une moyenne de 1 et une variance de 0,5. Ici, le coefficient de correction est

r  =\frac {q(θ^{t})f(a) }{q(a)f(θ^{t})}

Et. Le programme est ci-dessous (l'importation de la bibliothèque etc. n'est pas écrite. Veuillez consulter le lien texte intégral plus tard).

theta = []

# Initial value
current = 3
theta.append(current)

n_itr = 100000

for i in range(n_itr):
    #Génération aléatoire à partir de la distribution proposée
    a = rand_prop()
    if a < 0:
        theta.append(theta[-1])
        continue

    r = (q(current)*f_gamma(a)) / (q(a)*f_gamma(current))
    assert r > 0
    
    if r < 0:
        #reject
        theta.append(theta[-1])
        continue
    if r >= 1 or r > st.uniform.rvs():#Générez un nombre aléatoire et acceptez-le s'il est supérieur à r
        # Accept
        theta.append(a)
        current = a
    else:
        #Reject
        theta.append(theta[-1])

Le point est de savoir comment exprimer s'il faut accepter avec une probabilité de $ r $.

if r >= 1 or r > st.uniform.rvs()

Il est décidé que la méthode st.uniform.rvs () générera des nombres aléatoires dans une distribution uniforme de 0 $ à 1 $ et acceptera les plus grands.

À la suite de ce programme, vérifiez si la distribution postérieure peut être reproduite.


plt.title("Histgram from MCMC sampling.")
plt.hist(theta[n_burn_in:], bins=30, normed=True)

xx = np.linspace(0, 2.5,501)
plt.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))
plt.show()

017_100000(loc=1).png

Les données affichées en bleu sous forme d'histogramme sont la distribution de fréquence de $ θ $ déplacée par la méthode MH. Il est normalisé et corrigé pour que la somme soit 1. La couleur orange étant la distribution postérieure que nous voulions à l'origine reproduire, nous pouvons confirmer qu'elle peut être reproduite.

Eh bien, il semble que cela ait été fait sans aucun problème, mais dans la distribution proposée, la moyenne était de 1 et la variance de 0,5. Ceci est arbitraire et peut ne pas fonctionner. Par exemple, si une distribution normale avec une moyenne de 3 et une variance de 0,5 est la distribution proposée, 018.png

Ça va changer comme ça. En fait, dans la méthode MH indépendante, il y a un problème que la distribution ne converge pas bien à moins que la distribution proche de la distribution cible est une distribution stable. En d'autres termes, la méthode MH indépendante fait une grande différence dans les résultats jusqu'à convergence en fonction de la qualité de la distribution proposée. Cette fois, la distribution a posteriori a été clarifiée par souci de simplicité, mais dans de nombreux cas, la distribution est inconnue dans l'analyse réelle des données. Par conséquent, nous devons réfléchir à des moyens de résoudre ce problème.

Méthode MH de marche aléatoire

Un algorithme appelé méthode MH de marche aléatoire a été proposé comme méthode pour résoudre ce problème.

Candidats aux suggestions

a =θ^{t} + e

ça ira. Pour $ e $, une distribution normale avec une moyenne de 0 ou une distribution uniforme est utilisée.

Si vous choisissez une distribution symétrique telle qu'une distribution normale ou une distribution uniforme pour la distribution proposée,

q(a|θ^{t}) = q(θ^{t}|a)

Et la distribution proposée. Par conséquent, le coefficient de correction r dans la méthode de marche aléatoire MH est

r = \frac {f(a)}{f(θ^{t})}

Et la distribution proposée disparaît toujours.

Mettre en œuvre et vérifier (méthode MH de marche aléatoire)

Maintenant, implémentons et confirmons le même problème qu'avant.

# MCMC sampling

theta = []
current  = 5
theta.append(current)
n_itr = 100000
prop_m = 5
prop_sd = np.sqrt(0.10)
x = np.linspace(-1,5,501)

def rand_prop1(prop_m):
    return st.norm.rvs(prop_m, prop_sd)

for i in range(n_itr):
    a = rand_prop1(current)  #Génération aléatoire à partir de la distribution proposée
    r = f_gamma(a) / f_gamma(current)
    
    if a < 0 or r < 0:
        #reject
        theta.append(theta[-1])
        pass
    elif r >= 1 or r > st.uniform.rvs():
        # Accept
        theta.append(a)
        current = a
        status = "acc"
        col = "r"
    else:
        #Reject
        theta.append(theta[-1])
        pass

La description du graphique est ci-dessous.

plt.title("Histgram from MCMC sampling.")
n, b, p = plt.hist(theta[n_burn_in:], bins=50, density=True)

xx = np.linspace(0, 2.5,501)
plt.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))
plt.show()

kk = 11.
tt = 13.
print ("sample mean:{0:.5f}, sample std:{1:.5f}".format(np.mean(theta), np.std(theta)))
#Valeur théorique attendue / dispersion
print("mean:{0:.5f}, std:{1:.5f}".format( kk/tt, np.sqrt(kk*((1/tt)**2))) )

019.png

Il s'est avéré que la distribution était parfaitement reproduite. Comparé à la méthode MH indépendante que j'ai mentionnée plus tôt, j'ai trouvé que cet algorithme de méthode MH de marche aléatoire semble être sûr à utiliser.

Cependant, le problème de la méthode MH

Je voudrais dire que c'est un soulagement, mais en réalité, il y a encore des problèmes à résoudre. Après tout, vous devez spécifier correctement la valeur de $ e $, donc cela nécessite également un concept d'hyper-paramètre.   Si la variance est faible, le pourcentage de transitions acceptées est élevé, mais l'avance dans l'espace d'états est une marche aléatoire lente, qui nécessite un long temps de convergence. En revanche, plus la dispersion est élevée, plus le taux de rejet est élevé.   La méthode hamiltonienne de Monte Carlo est un moyen de résoudre ce problème. J'espère que cela sera résumé la prochaine fois.

À la fin

Cette fois, nous avons résumé la méthode Metropolis Hasting. Penser comment converger vers la moyenne de la distribution postérieure revient à optimiser les paramètres de poids dans un réseau neuronal.

Après tout, bien que les algorithmes soient différents, le but que je veux faire est similaire, alors j'ai pensé que si je pouvais comprendre cela, je pourrais voir le plaisir et l'ampleur des statistiques et de l'apprentissage automatique.

Le texte intégral du programme est stocké ici. https://github.com/Fumio-eisan/Montecarlo20200516

Recommended Posts

Comprendre la méthode Metropolitan Hasting (une des méthodes de la méthode Monte Carlo en chaîne de Markov) avec implémentation
La première méthode de Monte Carlo en chaîne de Markov par PyStan
Augmentez la vitesse de la méthode Monte Carlo de l'implémentation de découpage Cython.
[Statistiques] Visualisez et comprenez la méthode Hamiltonian Monte Carlo 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]
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
Essayez d'implémenter la méthode Monte Carlo en Python
Calcul de l'itinéraire le plus court selon la méthode de Monte Carlo
Sauvegardez la sortie de GAN une par une ~ Avec l'implémentation de GAN par PyTorch ~
Trouvez le ratio de la superficie du lac Biwa par la méthode de Monte Carlo
Simuler la méthode Monte Carlo en Python
4 méthodes pour compter le nombre d'occurrences d'entiers dans un certain intervalle (y compris la méthode imos) [implémentation Python]
Comprendre les images de diverses opérations matricielles utilisées dans Keras (Tensorflow) avec des exemples
En voici une, je vais résumer les applications équipées "d'intelligence artificielle" qui m'intéressaient