Блог

Я закодировал простой алгоритм 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)
  

Когда я наношу теты на гистограмму, я получаю следующую форму задней части … что определенно неверно.

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

Кажется, что выборки в хвостах более вероятны под задней частью. Но почему?

Я осознал свою ошибку. Я не должен был брать журнал коэффициентов ненормализованного заднего.