Markov Chain Monte Carlo and Metropolis Hasting

Markov Chain Monte Carlo (MCMC) is a sampling process where the next sample depends on the previous. Metropolis hasting is an accept-reject rule that indicates which random sample to keep or discard.

Markov Chain Monte Carlo

Markov Chain is a chain process that the next outcome is based on previous.

Monte Carlo is a random sampling process where repeatedly random sample to achieve a certain result. For example, if we repeatedly sample 10000 samples from a normal distribution — N (0,1). Then plot the histogram of 10000 samples, it will look like the initial distribution N(0,1).

import numpy as np
import matplotlib.pyplot as plt
# Monte Carlo
samples = np.random.normal(0, 1, 10000)
figure, ax = plt.subplots(1,2)
ax[0].hist(samples)
ax[0].set_title('Histogram of 10k samples from N(0,1)')
ax[1].plot(samples)
ax[1].set_title('10k Random samples')
plt.show()

MCMC is a special type of Monte Carlo, it is random sampling process but the next sample is depend on the previous (Markov Chain). For example, we repeatedly sample 10000 samples from N(0,1) but we replace the mean of the normal distribution with the previous sample:

Iteration 1: Sample x1 from N(0,1)
Iteration 2: Sample x2 from N(x1,1)
...
Iteration N: Sample xn from N(xn-1, 1)

The simulation in python:

import numpy as np
import matplotlib.pyplot as plt
# Markov Chain Monte Carlo
sample_t1 = np.random.normal(0, 1, 1)
samples = [sample_t1]
for t in range(1,10000):
    samples_t = np.random.normal(samples[t-1], 1, 1)
    samples = np.append(samples, samples_t)

figure, ax = plt.subplots(1,2)
ax[0].hist(samples)
ax[0].set_title('Histogram of 10k samples')
ax[1].plot(samples)
ax[1].set_title('10k Random samples')
plt.show()

This process is also called ‘random walk’, that the ending distribution is not always the same as the initial distribution. In this case, the histogram of 10k samples is not representing N(0,1).

Therefore, we might want to accept or reject some of the samples through the process to be able to get back to the initial distribution. Just like when take a random walk, we will need to accept or reject some routes to be able to get back to the starting point. Metropolis Hasting is one type of the accept or reject rule.

Metropolis Hasting

Let’s use the coin toss as an example.

I have a coin, and I want to test if it an unbiased coin. We have one parameter p(Head) to define. p(Head) is the probability of head. The test hypothesis is:

Hypothesis: 

H0: p(Head) = θ = 0.5
H1: p(Head) = θ is not 0.5

Option 1. T test: 
if p(t>t0.975) + p(t<t0.025) > 0.05, then accept null hypothesis,
else, accept althernative hypothesis

Option 2. Use quantile of θ:
if θ=0.5 is within 2.5% quatile and 97.5% quantile, then accept null hypothesis,
else, accept alternative hypothesis

There are three probabilities in Metropolis Hasting — Prior, Likelihood, Posterior.

Prior

Prior is the probability of p(Head) is a certain value, written as p(p(Head)=θ) or p(θ). It is determined by the prior distribution. For example, say this is a brand new coin and I don’t have any prior knowledge, I can use a uniform distribution — U(0,1) as prior distribution, since p(Head) can be any number in [0,1]. Randomly drew θ from U(0,1) and θ = 0.43, then the prior p(θ) is p(0.43) = 1 since f(x) =1 for all x in [0,1] for Uniform(0,1).

Note U(0,1) equals Beta(1,1).

Prior distribution

Likelihood

Likelihood is the probability of getting certain data when p(Head) =θ, written as p(data|θ). It is determined by the distribution of the data. In this case, I repeatedly toss my coin for 10 times, and record head and tail each toss. Each toss is Bernoulli trial with two outcomes Head or Tail. The data is a binomial distribution. The likelihood is p(data|θ) = Binomial(k, n, θ). k is number of heads, n is total toss (n=10 in this example), θ is p(Head).

Posterior

Posterior is the probability of getting p(Head) = θ given data, written as: p(θ|data). It can be calculated using prior and likelihood.

p(data) is not associated with parameter-θ, we can think it as a constant, and simplify this equation into:

Metropolis hasting

The accept and reject rule of Metropolis hasting is:

Coin Toss Experiment

The algorithm is:

# initial
1. Drew a θ from uniform(0,1)
       For example: θ = 0.07
2. Compute prior 
       p(θ=0.07) = 1
3. Toss the coin for 10 times: H-heads, T-Tails and count heads 
       For example: [H, H, T, T, T, H, H, T, T, T]
       Total number of Heads = 4
4. Compute likelihood:
       Likelihood = Binomial (4, 10, 0.07) = 0.003
5. Compute posterior:
       posterior = prior * likelihood = 1 * 0.003 = 0.003

# iteration 1
6. Drew a θ from uniform(0,1)
       For example: θ = 0.43
7. Compute prior 
       p(θ=0.43) = 1
8. Toss the coin for 10 times: H-heads, T-Tails and count heads 
       For example: [H, H, H, T, T, H, H, T, T, T]
       Total number of Heads = 5
9. Compute likelihood:
       Likelihood = Binomial (5, 10, 0.43) = 0.229
10. Compute posterior:
       posterior = prior * likelihood = 1 * 0.229 = 0.229
11. Compute ratio of posterior(θt) / posterior(θt-1)
       r = 0.229/0.003 = 7.6
12. Accept or reject:
       r > 1, accept θ = 0.43, update prior =1, 
                               update likelihood = 0.229, 
                               update posterior = 0.229

In the following table, I wrote down the first 7 iterations. Data is the total head counts of 10 toss. u is the random sample from uniform distribution (0,1) only when r <1.

From the first 7 iteration, I collected few θ: 0.43, 0.54, 0.62, 0.32.

Get Keru Chen’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

I will do this experiment N times until I have 1000 of θ. N≥1000 since I am rejecting θ along the way.

# MH
import scipy
import random
import numpy as np
#  initial
theta0 = np.random.beta(1,1,1)
prior0 = scipy.stats.beta.pdf(theta0, 1, 1)
head_counts_0 = sum(np.random.randint(2, size=10))
likelihood_0 = scipy.stats.binom.pmf(head_counts_0, 10, theta0)

# collect
theta_accept = theta0
prior_accept = prior0
likelihood_accept = likelihood_0

i = 0
while len(theta_accept) <= 1000:
    theta = np.random.beta(1,1,1)
    prior = scipy.stats.beta.pdf(theta, 1, 1)
    head_counts = sum(np.random.randint(2, size=10))
    likelihood = scipy.stats.binom.pmf(head_counts, 10, theta)
    ratio = (prior*likelihood) / (prior_accept[i] * likelihood_accept[i])
    if ratio >=1:
        theta_accept = np.append(theta_accept, theta)
        prior_accept = np.append(prior_accept, prior)
        likelihood_accept = np.append(likelihood_accept, likelihood)
    else:
        u = np.random.beta(1,1,1)
        if ratio >= u:
            theta_accept = np.append(theta_accept, theta)
            prior_accept = np.append(prior_accept, prior)
            likelihood_accept = np.append(likelihood_accept, likelihood)
        else:
            print('ratio is: ', ratio)
            print('u is: ', u)

# drop the first theta, prior, likelihood
theta_accept = theta_accept[1:]
prior_accept = prior_accept[1:]
likelihood_accept = likelihood_accept[1:]

Here is the first 20 θ that is accepted:

theta_accept[:20]
Out[6]: 
array([0.4182804 , 0.27898128, 0.87778638, 0.56595161, 0.70156645,
       0.49273126, 0.47311259, 0.75905831, 0.40404581, 0.4929927 ,
       0.43557259, 0.8074565 , 0.50941352, 0.26218544, 0.14406795,
       0.32712108, 0.6867593 , 0.54402593, 0.52818101, 0.4457688 ])

Plot the histogram of θ and this will be posterior distribution:

plt.hist(theta_accept)
plt.show()

Posterior distribution of θ given data

Now going back to our hypothesis, set significant level ɑ = 0.05, and this is a two-side test:

H0: p(Head) = θ = 0.5, this is an unbiased coin
H1: p(Head) = θ is not 0.5, this is a biased coin

Option 1. T test: 
if p(t>t0.975) + p(t<t0.025) > 0.05, then accept null hypothesis;
else, accept althernative hypothesis

Option 2. Use quantile of θ:
if θ=0.5 is in 2.5% quatile and 97.5% quantile, then accept null hypothesis;
else, accept alternative hypothesis

Note that the degree of freedom (df) is 10, since we toss 10 times each iteration. Given data:

# t test
mean = np.mean(theta_accept)
print('mean of theta, ', mean)
sd = np.std(theta_accept)
print('std of theta, ', sd)
df = 10
se = sd / np.sqrt(df)
print ('standard error of theta, ', se)
t = np.abs((mean-0.5) / se)
print('t = ', t)

p = (1-scipy.stats.t.cdf(t, 10))*2
print('p value from t test is, ', p)

# output 
mean of theta,  0.4967416562316055
std of theta,  0.2359528875945374
standard error of theta,  0.07461485452924263
t =  0.04366883496526178
p value from t test is,  0.9660280570824347

P value is greater than 0.05, fail to reject null hypothesis: p(Head) = θ = 0.5, this is an unbiased coin.

The second option is using quantile. First sort the 1000 θ, the 2.5% quantile will be 25th θ after sorting and 97.5% quantile will be 975th θ after sorting.

# sort theta
theta_accept_sorted = np.sort(theta_accept)
theta_25_quantile = theta_accept_sorted[25]
theta_975_quantile = theta_accept_sorted[975]

print('the 2.5% quantile is: ', theta_25_quantile)
print('the 97.5% quantile is: ', theta_975_quantile)

# output
the 2.5% quantile is:  0.08962180547280807
the 97.5% quantile is:  0.9119102420388027

Therefore, 0.5 is within [0.089, 0.912], fail to reject null hypothesis: p(Head) = θ = 0.5, this is an unbiased coin.

It is a true unbiased coin because I used “np.random.randint(2, size=10)” function to generate head and tail. This function uses a discrete uniform distribution U[0,1], hence the probability for p(x=0)=p(x=1)=0.5.

Last let’s visualize the prior distribution — initial distribution of θ, the likelihood distribution, and posterior distribution —final distribution of θ given data by plotting histogram of samples.

# plot
theta_prior = np.random.beta(1,1,1000)
likelihood = likelihood_accept
theta_posterior = theta_accept

plt.hist(theta_prior, histtype='step', edgecolor='blue', linewidth=2)
plt.hist(likelihood, histtype='step', edgecolor='orange', linewidth=2)
plt.hist(theta_posterior, histtype='step', edgecolor='green', linewidth=2)
plt.legend(['prior distribution', 'likelihood distribution', 'posterior distribution'])
plt.show()

From histogram we can see that the prior distribution is a uniform distribution from [0,1] or beta(1,1). And posterior distribution is drifted by the data (represented by the orange line) to a different Beta (a,b) distribution. Beta prior is a conjugate prior, it shares the same distribution with the posterior. And the posterior distribution can be obtain through prior distribution and likelihood distribution. In this case, the prior distribution is Beta(1,1) and the likelihood distribution is Binomial(k, n, θ), then we can obtain the posterior distribution is Beta(1+k, 1+(n-k)).

That’s all for the MCMC and MH :)

Hope you enjoyed reading :)

Reference

Special thanks to Youtube — StataCorp LLC. I have learned a lot from their videos about Bayesian statistics.