Metropolis hastings算法用于简单抛硬币

Question

我已经编写了一个简单的MH算法来估计抛硬币落地头部的后验概率。我已经看过代码很多次了，但我不明白大都会算法为什么不能正常工作。

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, trials, heads, prior):
   log_likelihood_theta_curr = log_likelihood(theta_curr, trials, heads)
    log_likelihood_theta_prop = log_likelihood(theta_prop, trials, 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, trials, heads, 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, theta_prop, trials, heads, 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)

当我在直方图上绘制thetas时，我得到了如下的后部形状...这肯定是不正确的。

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

​

​

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

Answer 1

艾维意识到我的错误了。我不应该取非正规化后验概率的对数。