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 theta_curr_2 = theta_curr theta_prop_1 = theta_prop theta_prop_2 = theta_prop # 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.
$$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.