用于简单硬币翻转的Metropolis Hastings算法

问题描述

我已经编码了一个简单的MH算法，以估计抛硬币落下的可能性。我已经看过很多遍代码了，但我不明白Metropolis算法如何无法正常工作。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
import math

def nCr(n,r):
    f = math.factorial
    return f(n) / f(r) / f(n-r)

def log_likelihood(theta,trials,heads):
    likelihood = nCr(trials,heads) * (theta**heads) * (1 - theta) ** (trials - heads)
    return np.log(likelihood)

def log_prior(prior,theta):
    return np.log(prior.pdf(theta))

def compute_acceptance_prob(theta_curr,theta_prop,heads,prior):
   log_likelihood_theta_curr = log_likelihood(theta_curr,heads)
    log_likelihood_theta_prop = log_likelihood(theta_prop,heads)
    
    log_prior_curr = log_prior(prior,theta_curr)
    log_prior_prop = log_prior(prior,theta_prop)
    
    log_post_curr = log_likelihood_theta_curr + log_prior_curr
    log_post_prop = log_likelihood_theta_prop + log_prior_prop
    
    
    return min(1,log_post_prop / log_post_curr)



def Metropolis_hastings(prior,a,b,MAX_ITERS = 100000):
    thetas = []
    theta_curr = beta.rvs(a,b)
    for i in range(MAX_ITERS):
        theta_prop = theta_curr + np.random.normal(0,0.05)
        #print(theta_prop)
        if theta_prop > 1 or theta_prop < 0:
            theta_prop = theta_curr
        
        prob = compute_acceptance_prob(theta_curr,prior)
        r = np.random.uniform(0,1)
        
        if prob > r:
            theta_curr = theta_prop
        
        thetas.append(theta_curr)
    return thetas

thetas = Metropolis_hastings(theta_1,5,6,20,10)

当我在直方图上绘制theta时，我得到了以下后验形状...这绝对是不正确的。

plt.hist(thetas[20000:],histtype='stepfilled',color = 'darkred',bins=30,alpha=0.8,density=True);

似乎在尾部下方的尾部样本更有可能。但是为什么呢？

解决方法

我意识到我的错误。我不应该记录未归一化后验的比率。