# Metropolis-Hastings Sampler

## Goals

In this tutorial, we look at the Metropolis-Hastings sampler. This tutorial will cover how to:

1. Implement the Metropolis-Hastings sampler using two different variations of acceptance probabilities:
• An independence chain
• A random walk chain
2. Calculate the distribution's mean and variance-covariance matrix.

## Introduction

The Metroplis-Hastings sampler is an iterative algorithm which produces a sequence of draws $\theta^{(r)}$ which can be used to estimate sample parameters such that

$$\hat{g} = \frac{\sum_{r=1}^{R} g(\theta^{(r)})}{R}$$

The intuition is relatively simple. The algorithm:

• Picks a candidate draw from the specified candidate generating density.
• Determines if the draw is accepted or rejected based on a computed acceptance probability.
• For each draw $r = 1, ..., R$ sets $\theta^{(r)}$ equal to $\theta^{(r-1)}$, if the draw is rejected.

In this particular tutorial, we will examine a single parameter $\theta$ with the posterior distribution given by

$$p(\theta|y) \propto exp\bigg(-\frac{1}{2}|(\theta)|\bigg)$$

### The Metropolis-Hastings steps

The general Metropolis-Hastings algorithm can be broken down into simple steps:

1. Set up sampler specifications, including number of iterations and number of burn-ins draws.
2. Choose a starting value for $\theta$.
3. Draw $\theta^*$ from the candidate generating density.
4. Calculate the acceptance probability $\alpha(\theta^{(r-1)}, \theta^*)$.
5. Set $\theta^{(r)} = \theta^*$ with probability $\alpha(\theta^{(r-1)}, \theta^*)$, or else set $\theta^{(r)} = \theta^{(r-1)}$.
6. Repeat steps 3-5 for $r = 1, ... , R$.
7. Average the draws to produce the Metropolis-Hastings estimates for posterior functions of interest (mean, variance, etc).

## The independent chain

### The candidate generating density

We will first consider an independent chain which uses a normal distribution to generate the candidate draws

$$N\big(0,\text{ }d^2\big)$$

### The acceptance probability

The acceptance probability for this algorithm depends on the previous draw, $\theta^{(r-1)}$, the candidate draw, $\theta^*$, and the performance parameter, $d$ :

$$\alpha\big(\theta^{(r-1)}, \theta^*\big) = min \bigg[exp\Big[\frac{1}{2}\Big( | \theta^{(r-1)}| - | \theta^*| - \bigg(\frac{\theta^{(r-1)}}{d}\bigg)^2 + \bigg(\frac{\theta^*}{d}\bigg)^2\Big) \Big],1\bigg]$$

### Metropolis-Hastings settings

Our first step is to set up the specifications of the Metropolis-Hastings algorithm. We will :

• Keep 10,000 total draws
• Set a burn-in sample of 100 draws.
• Set the starting values of $\theta$ to 1.
• Set the standard deviation of the candidate draw to 6.
new;

/*
** Set up Metropolis Hastings chain
*/
burn_in = 100;

// Specify the total number of replications
keep = 10000;

/*
** Standard deviation for the increment
** in the independent chain M-H algorithm
*/
sd_ic = 6;

// Initialize MH sums to zero
theta_mean_ic = 0;
th2_mo_ic = 0;
theta_draw_ic = zeros(keep+burn_in, 1);

// Specify starting value for independent chain theta draw
theta_draw_ic[1, 1] = 1;

// Set count of the number of accepted draws
count_ic = 0;

### Implementing the algorithm

Next, we will implement the Metropolis-Hastings algorithm using a for loop. Within each iteration of our loop we need to :

1. Draw a candidate $\theta$ using the GAUSS rndn function.
2. Compute the acceptance probability.
3. Set $\theta^{(r)} = \theta^*$ with probability $\alpha(\theta^{(r-1)}, \theta^*)$, or else set $\theta^{(r)} = \theta^{(r-1)}$.
// Set random seed for repeatable random numbers
rndseed 97980;

for i(2, burn_in+keep, 1);
// Front multiply by sd_ic
theta_can_ic = sd_ic*rndn(1, 1);

// Calculate acceptance probability
acc_prob_ic = minc(exp(.5*(abs(theta_draw_ic[i-1, 1]) - abs(theta_can_ic) + (theta_can_ic/sd_ic)^2 - (theta_draw_ic[i-1, 1]/sd_ic)^2))|1);

// Accept candidate draw with probability acc_prob_ic
if acc_prob_ic > rndu(1, 1);
theta_draw_ic[i, 1] = theta_can_ic;
count_ic = count_ic + 1;
else;
theta_draw_ic[i, 1] = theta_draw_ic[i-1, 1];
endif;

/*
** Store remainder for computing
** sample parameters.
*/
if i>burn_in;
/*
** This keeps a running sum of
** of theta. It can be used to
** find the mean of theta.
*/
theta_mean_ic = theta_mean_ic + theta_draw_ic[i, 1];
/*
** This keeps a running sum of
** of square of theta. It will
** be used to find the standard
** deviation of theta.
*/
th2_mo_ic = th2_mo_ic + theta_draw_ic[i, 1].^2;
endif;
endfor;

## The random walk chain

As an alternative to the independent chain, let's also consider the random walk chain. The random walk chain Metropolis-Hastings algorithm generates candidate draws from the process

$$\theta^* = \theta^{(r-1)} + z$$

where $z$ is a random variable that follows a normal increment.

### The candidate generating density

The iteration candidate, $\theta^*$ is drawn from $N(\theta^{(r-1)},\text{ } c^2)$ where, $c$ is a performance parameter. Note how this differs from the independent chain which is drawn from $N(0,\text{ } c^2)$.

### The acceptance probability

In the random walk case, the acceptance probability is given by

$$\alpha\big(\theta^{(r-1)}, \theta^*\big) = min \bigg[exp\Big[\frac{1}{2}\Big( | \theta^{(r-1)}| - | \theta^*|\Big) \Big],1\bigg]$$

### Implementing the algorithm

We will again use a for loop to complete the steps of our Metropolis-Hastings sampler

// Set random seed
rndseed 10385;

// Standard deviation for RW M-H algorithm
sd_rw = 4;

// Initialize MH variables
theta_mean_rw = 0;
th2_mo_rw = 0;
theta_draw_rw = zeros(keep+burn_in, 1);

// Specify starting value for theta draw
theta_draw_rw[1, 1] = 1;

// Set count of the number of accepted draws
count_rw = 0;

for i(2, burn_in+keep, 1);
// Draw a candidate for random walk chain
theta_can_rw = theta_draw_rw[i, 1] + sd_rw*rndn(1, 1);

// Calculate acceptance probability
acc_prob_rw = minc(exp(.5*(abs(theta_draw_rw[i, 1]) - abs(theta_can_rw)))|1);

// Accept candidate draw with probability acc_prob_rw
if acc_prob_rw > rndu(1, 1);
theta_draw_rw[i, 1] = theta_can_rw;
count_rw = count_rw+1;
else;
theta_draw_rw[i, 1] = theta_draw_rw[i-1, 1];
endif;

/*
** Store remainder for computing
** sample parameters.
*/
if i>burn_in;
/*
** This keeps a running sum of
** of theta. It can be used to
** find the mean of theta.
*/
theta_mean_rw = theta_mean_rw + theta_draw_rw[i,1];
/*
** This keeps a running sum of
** of square of theta. It will
** be used to find the standard
** deviation of theta.
*/
th2_mo_rw = th2_mo_rw + theta_draw_rw[i, 1].^2;
endif;
endfor;

## Compute sample statistics

Finally, we use our stored data to find our sampler mean and variance.

print "Number of burn in replications";
burn_in;

print "Number of included replications";
keep;

print "Proportion of accepted candidate draws: Independence chain M-H";
count_ic/(keep+burn_in);

print "Proportion of accepted candidate draws: Random walk chain M-H";
count_rw/(keep+burn_in);

theta_mean = (theta_mean_ic~theta_mean_rw)./keep;
th2_mo = (th2_mo_ic~th2_mo_rw)/keep;
thvar = th2_mo - theta_mean.^2;

print "Posterior Mean and Variance, Ind Chain then RW Chain";
theta_mean;

thvar;

The above code should print the following output:

Number of burn in replications
100.0000
Number of included replications
10000.0000
Proportion of accepted candidate draws: Independence chain M-H
0.48455446
Proportion of accepted candidate draws: Random walk chain M-H
0.52603960
Posterior Mean and Variance, Ind Chain then RW Chain
0.038943707      0.075589115
8.1069691        8.0770516 

### Conclusion

Congratulations! You have:

• Implemented the Metropolis-Hastings sampler using as an independent chain.
• Implemented the Metropolis-Hastings sampler using as a random walk chain.
• Calculated the distribution's mean and variance-covariance matrix.

For your convenience, the entire code is below.

new;

/*
** Set up Metropolis Hastings chain
*/
burn_in = 100;

// Specify the total number of replications
keep = 10000;

/*
** Standard deviation for the increment
** in the independent chain M-H algorithm
*/
sd_ic = 6;

// Initialize MH sums to zero
theta_mean_ic = 0;
th2_mo_ic = 0;
theta_draw_ic = zeros(keep+burn_in, 1);

/*
** Specify starting value for
** independent chain theta draw
*/
theta_draw_ic[1, 1] = 1;

// Set count of the number of accepted draws
count_ic = 0;

// Set random seed for repeatable random numbers
rndseed 97980;

for i(2, burn_in+keep, 1);
// Candidate draw from normal distribution
theta_can_ic = sd_ic*rndn(1, 1);

// Calculate acceptance probability
acc_prob_ic = minc(exp(.5*(abs(theta_draw_ic[i-1, 1]) - abs(theta_can_ic) +
(theta_can_ic/sd_ic)^2 - (theta_draw_ic[i-1, 1]/sd_ic)^2))|1);

/*
** Accept candidate draw with probability
** theta_can_ic
*/
if acc_prob_ic > rndu(1, 1);
theta_draw_ic[i, 1] = theta_can_ic;
count_ic = count_ic + 1;
else;
theta_draw_ic[i, 1] = theta_draw_ic[i-1, 1];
endif;

if i>burn_in;
/*
** This keeps a running sum of
** of theta. It can be used to
** find the mean of theta.
*/
theta_mean_ic = theta_mean_ic + theta_draw_ic[i, 1];

/*
** This keeps a running sum of
** of square of theta. It will
** be used to find the standard
** deviation of theta.
*/
th2_mo_ic = th2_mo_ic + theta_draw_ic[i, 1].^2;
endif;
endfor;

// Set random seed for repeatable random numbers
rndseed 10385;

/*
** Standard deviation for the
** RW M-H algorithm
*/
sd_rw = 4;

// Initialize MH sums to zero
theta_mean_rw = 0;
th2_mo_rw = 0;
theta_draw_rw = zeros(keep+burn_in, 1);

// Specify starting value for random walk chain theta draw
theta_draw_rw[1, 1] = 1;

// Set count of the number of accepted draws
count_rw = 0;

for i(2, burn_in+keep, 1);
// Draw a candidate for random walk chain
theta_can_rw = theta_draw_rw[i-1, 1] + sd_rw*rndn(1, 1);

// Calculate acceptance probability
acc_prob_rw = minc(exp(.5*(abs(theta_draw_rw[i-1, 1]) - abs(theta_can_rw)))|1);

// Accept candidate draw with probability acc_prob_rw
if acc_prob_rw > rndu(1, 1);
theta_draw_rw[i, 1] = theta_can_rw;
count_rw = count_rw+1;
else;
theta_draw_rw[i, 1] = theta_draw_rw[i-1, 1];
endif;

if i>burn_in;
/*
** This keeps a running sum of
** of theta. It can be used to
** find the mean of theta.
*/
theta_mean_rw = theta_mean_rw + theta_draw_rw[i,1];
/*
** This keeps a running sum of
** of square of theta. It will
** be used to find the standard
** deviation of theta.
*/
th2_mo_rw = th2_mo_rw + theta_draw_rw[i, 1].^2;
endif;
endfor;

print "Number of burn in replications";
burn_in;

print "Number of included replications";
keep;

print "Proportion of accepted candidate draws: Independence chain M-H";
count_ic/(keep+burn_in);

print "Proportion of accepted candidate draws: Random walk chain M-H";
count_rw/(keep+burn_in);

// Calculate mean using stored sum of theta
theta_mean = (theta_mean_ic~theta_mean_rw)./keep;

// Calculate standard deviation
th2_mo = (th2_mo_ic~th2_mo_rw)/keep;
thvar = th2_mo - theta_mean.^2;

print "Posterior Mean and Variance, Ind Chain then RW Chain";
theta_mean;

thvar;

### Have a Specific Question?

Get a real answer from a real person

### Need Support?

Get help from our friendly experts.