用于简单硬币翻转的 Metropolis 黑斯廷斯算法

Question

我编写了一个简单的 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, 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

我意识到我的错误。我不应该记录非标准化后验的比率。