Metropolis Hastings algorithm for joint posterior of probability of heads for 2 coins

Question

I am trying to implement a simple metropolis hastings algorithm to simulate the joint posterior of the probability of flipping heads for 2 coins, $\theta_1,\theta_2$. I am following the problem statement outlined in this tutorial by Patrick Ford.

I am not able to get why the algorithm keeps getting values that exceed the range of 0 < $\theta$ < 1. The metropolis Hastings algorithm should reject samples that are far away from the true posterior probability and accept samples that are closer to the true posterior.

This is the code that I have for the simple problem. Am I calculating the unnormalized posterior wrongly or something ?

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


# parameters, lets say theta_1 is unbiased and prior to our data we observed 5 heads and 6 tails.
theta_1 = beta(5,6)

# lets say theta_2 is biased. Prior to our data we observed 10 heads and 4 tails
theta_2 = beta(10,4)

# for theta_1, the data we obtain is 10 heads and 10 tails 
# for theta_2, the data we obtain is 15 heads and 6 tails

# start metropolis hastings. Posterior can be computed analytically as well.

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

def log_likeilihood(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):
    # likelihood * prior. Compute the log unnormalized posterior
    # logPost = logLLH + logPrior
    theta_curr_1 = theta_curr[0]
    theta_curr_2 = theta_curr[1]
    
    theta_prop_1 = theta_prop[0]
    theta_prop_2 = theta_prop[1]
    
#     print(theta_prop_1)
#     print(theta_prop_2)
#     print(theta_curr_1)
#     print(theta_curr_2)
    
    log_likelihood_theta_curr_1 = log_likeilihood(theta_curr_1, 20, 10) # 20 trials, 10 heads
    log_likelihood_theta_curr_2 = log_likeilihood(theta_curr_2, 21, 15) # 21 trials, 15 heads
    
    log_likelihood_theta_prop_1 = log_likeilihood(theta_prop_1, 20, 10) # 20 trials, 10 heads
    log_likelihood_theta_prop_2 = log_likeilihood(theta_prop_2, 21, 15) # 21 trials, 15 heads
    
    log_prior_curr_1 = log_prior(theta_1, theta_curr_1)
    log_prior_curr_2 = log_prior(theta_2, theta_curr_2)
    
    log_prior_prop_1 = log_prior(theta_1, theta_prop_1)
    log_prior_prop_2 = log_prior(theta_2, theta_prop_2)
    
    joint_likelihood_curr = log_likelihood_theta_curr_1 + log_likelihood_theta_curr_2
    joint_likelihood_prop = log_likelihood_theta_prop_1 + log_likelihood_theta_prop_2
    
    log_post_curr = joint_likelihood_curr + log_prior_curr_1 + log_prior_curr_2
    log_post_prop = joint_likelihood_prop + log_prior_prop_1 + log_prior_prop_2
    
    return min(1,log_post_prop / log_post_curr)

theta_curr = np.array([np.random.beta(5,6), np.random.beta(10,4)]) # sample from beta 
print(theta_curr)
MAX_ITERS = 10000
prop_mean = [0,0]
prop_cov = [[0.005,0],[0,0.005]]
thetas = []
for i in range(MAX_ITERS):
    # sample from gaussian proposal 
    #print(np.random.multivariate_normal(prop_mean, prop_cov))
    theta_prop = theta_curr + np.random.multivariate_normal(prop_mean, prop_cov)
    print(theta_prop)
    # evaluate unnormalized posterior
    prob = compute_acceptance_prob(theta_curr, theta_prop)
    r = np.random.uniform(0,1)
    # evaluate unnormalized posterior
    thetas.append(theta_curr)
    if prob > r:
        theta_curr = theta_prop

Have been thinking about it for awhile now. Alot of the samples from the posterior are $< 0 $ and $>1$, which Im not sure why.

EDIT

$$p(D_1,D_2|\theta_1,\theta_2) = p(D_1|\theta_1)*p(D_2|\theta_2)$$ Therefore the joint log likelihood is $log\, p(D_1|\theta_1) + log \,p(D_2|\theta_2)$ as what I have computed.

EDIT I realised that there is a mistake in computing the acceptance ratio. I should not have used the unnormalised log ratio in the computation.

Answer 1

The Binomial likelihood function$$\ell(\theta)={n \choose x}\theta^x(1-\theta)^{n-x}$$is only defined when $0\le\theta\le 1$. Computing $\ell(1.4)$ or $\ell(-2.3)$ makes no sense as there is no Binomial experiment with a probability of $1.4$ or $-2.3$. For a parameter value outside the possible parameter range $(0,1)$, the likelihood need be set to zero.

The referenced tutorial has this safety step:

# If any of the values in thetas_prop is not in [0, 1.0] 
# then replace their value with
# the same element in thetas_curr, otherwise keep their values
thetas_prop <-mapply(FUN =function(x, y)ifelse((x>1|x<0), y, x),
    x = thetas_prop,y = thetas_curr)

meaning that the proposed values for $\theta_1$ and $\theta_2$ are always in the range $(0,1)$. The truncation in the proposal means that the Metropolis-Hastings acceptance ratio as given in the tutorial

# Accept/reject logic
p_accept_theta_prop =pmin(posterior_prop/posterior_curr, 1.0)

should be checked. Indeed, the proposal density is $$f(\xi|\theta)=\Phi((0,1)|\theta)\frac{\varphi(\xi-\theta)\mathbb I_{(0,1)}(\xi)}{\Phi((0,1)|\theta)}+\left[1-\Phi((0,1)|\theta)\right]\mathbb I_{\xi=\theta}$$ where $\varphi$ denotes the Normal density and $\Phi$ the Normal cdf. The acceptance ratio is then $$\dfrac{\pi(\xi|x)}{\pi(\theta|x)}\,\dfrac{f(\theta|\xi)}{f(\xi|\theta)}$$ and symmetry in the proposal leads to justify the ratio as given.

But a perfectly valid and straightforward alternative is to simply set the likelihood to zero when the proposed value is outside the range $(0,1)$. This clearly prevents the Markov chain to move there.

Note on the minor side that this example is somewhat absurd as the two components of $\theta$ are independent a posteriori and the simulated bivariate Markov chain is less efficient than simulated two independent Markov chains.