[PYTHON] 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.
-
[Statistiques de Bayes à partir des bases: une introduction pratique à la méthode hamiltonienne de Monte Carlo](https://www.amazon.co.jp/%E5%9F%BA%E7%A4%8E%E3%81%8B%E3 % 82% 89% E3% 81% AE% E3% 83% 99% E3% 82% A4% E3% 82% BA% E7% B5% B1% E8% A8% 88% E5% AD% A6-% E3% 83% 8F% E3% 83% 9F% E3% 83% AB% E3% 83% 88% E3% 83% 8B% E3% 82% A2% E3% 83% B3% E3% 83% A2% E3% 83% B3% E3% 83% 86% E3% 82% AB% E3% 83% AB% E3% 83% AD% E6% B3% 95% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E5% AE% 9F% E8% B7% B5% E7% 9A% 84% E5% 85% A5% E9% 96% 80-% E8% B1% 8A% E7% 94% B0-% E7% A7% 80% E6% A8% B9 / dp / 4254122128)
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é **.
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.
- ** Méthode de hachage Metropolis (expliquée cette fois) **
- Méthode de Monte Carlo hamiltonienne
- Donne un échantillon
- Échantillonnage de tranches
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,
Donc ce noyau de transition
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 deux coefficients de correction [tex: c, c '] et il est redondant
- Il est difficile de prendre des mesures probabilistes car le coefficient de correction n'est pas compris entre 0 ou plus et 1 ou moins.
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 devient
Algorithme MH
Maintenant, l'algorithme MH est ci-dessous.
- Utilisez la distribution proposée $ q (| θ ^ {t}) $ pour générer un nombre aléatoire $ a $.
- Déterminez les propositions suivantes.
q(a|θ^{t})f(θ^{t}) > q(θ^{t}|a)f(a)
- Vrai: À partir de la condition d'inégalité de l'équation 1, il est déterminé que $ θ = a $, $ θ '= θ ^ {t} $. Dans ce cas, la distribution proposée est $ q (θ | θ ') $, il est donc nécessaire de faire une correction probabiliste en utilisant le coefficient $ r $ sur le côté gauche de l'équation 2.
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 $.
- Faux: À partir de la condition d'inégalité de l'équation 1, il est déterminé que $ θ '= a $, $ θ = θ ^ {t} $. Dans ce cas, la distribution proposée est $ q (θ '| θ) $, il est donc nécessaire de faire une correction probabiliste en utilisant le coefficient $ r' $ sur le côté droit de l'équation 2. Acceptez toujours $ a $ avec une probabilité de $ 1 (= r ') $ et définissez $ θ ^ {t + 1} = a $.
- 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 (θ) $.
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()
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,
Ç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))) )
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