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

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)

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.

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.

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.

• Thank you for your help. I thought the chain should favour likelihood’s that are closer to the posterior distribution and hence theta values that are within 0-1 should be favoured over theta values that are >1 or <0? Or is there some undefined behaviour when likelihood’s > 1 and <0 ? Sep 12, 2020 at 8:06
• Also what do you mean by independent a posteriori ? Sep 12, 2020 at 8:07
• See update in the answer. Note that I exclude parameter values outside $(0,1)$, not likelihood values outside $(0,1)$. Sep 12, 2020 at 8:18
• When I run the code I get samples that are >1 and <0. Intuitively however shouldn’t those values be more likely to be rejected as they are far away from the true posterior, which lies between 0 - 1? Sep 12, 2020 at 8:24
• We have a communication problem: I am telling you values outside $(0,1)$ should be rejected 100% of the time, as being impossible, by extending the likelihood function or the prior density to be zero outside $(0,1)$. Sep 12, 2020 at 9:00