Metropolis-Hastings Sampler

Introduction

In this tutorial, we will

  1. Implement the Metropolis-Hastings sampler using two different variations of acceptance probabilities:
    • An independence chain
    • A random walk chain
  2. Compare the samples from the two methods in side-by-side histograms.
  3. Calculate the distribution's mean and variance-covariance matrix.

Background

The Metroplis-Hastings sampler is an iterative algorithm which produces a sequence of draws $\theta^{(r)}$ for $r = 1, ... , R$. The draws are averaged to calculate a parameter estimate such that

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

The algorithm picks a candidate draw from the specified candidate generating density. An acceptance probability then determines whether the candidate draw is accepted or rejected. If the candidate draw is rejected $\theta^{(r)}$ is set to $\theta^{(r-1)}$.

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 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 $p(\theta_1|y,\theta_2^{(0)})$.
  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).

We will first consider an independence chain with the $N\big(0,\text{ }d^2\big)$ candidate generating density. The acceptance probability has the form

$$ \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]$$

where $d$ is a performance algorithm. For this example we will set $d = 6$. However, in practice convergence diagnostics should be used to help choose $d$ to optimize the sampler.

Our first step is to set up the specifications of the Metropolis-Hastings algorithm. We will keep 10,000 total draws and set a burn-in sample of 100 draws. The starting values of $\theta$ are set to 1 and the standard deviation of the candidate draw is set to 6.

Metropolis-Hastings settings

new;

//Set up Metropolis Hastings chain
//Discard r0 burnin replications and 
burn_in = 100;

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

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

//Initialize MH sums to zero - we use a 2x1  
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;

Compute Independent Chain

Next, we will set up the independent chain for loop_. The for loop offers an efficient and compact method for algorithms that rely on loop for a given number of iterations. The GAUSS for keyword requires three inputs: start, stop, and step

//This illustrates general usage of a 'for' loop
//Do not run this code
for i(start, stop, step);
    //loop body commands
endfor;

Within our for 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);
    //Draw a candidate for independence chain
    //Front multiply by sd 
    theta_can_ic = sd_ic*rndn(1, 1);

    //Calculate acceptance probability
    acc_prob_ic = minc(exp(.5*(abs(theta_draw_ic[i, 1]) - abs(theta_can_ic) + (theta_can_ic/sd_ic)^2 - (theta_draw_ic[i, 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;
        //Discard r0 burn in draws
        theta_mean_ic = theta_mean_ic + theta_draw_ic[i, 1];  
        th2_mo_ic = th2_mo_ic + theta_draw_ic[i, 1].^2;
    endif;
endfor;

Compute Random Walk

As an alternative to the independent chain, we also will consider the random walk chain. In practice, it would be better to program the random and independent chain sampler in the same for loop to minimize the amount of looping in our program. However, for the sake of demonstration, the complete random walk sampler is below in its own for loop.

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 iteration candidate, $\theta^*$ is drawn from $N(\theta^{(r-1)},c^2)$

where, $c$ is a performance parameter. In this case we set $c = 4$.

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

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

//c parameter
//Standard deviation for the increment in the RW M-H algorithm
sd_rw = 4;

//Initialize MH sums to zero - we use a 2x1  
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] + 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 accrw
    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;
        //Discard r0 burnin draws
        theta_mean_rw = theta_mean_rw + theta_draw_rw[i,1];  
        th2_mo_rw = th2_mo_rw + theta_draw_rw[i, 1].^2;
    endif;
endfor;

Finally, we 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.2676
Proportion of accepted candidate draws: Random walk chain M-H
  0.3361
Posterior Mean and Variance, Ind Chain then RW Chain
  0.0337   0.0775
  8.0990   4.1984 

Plot results

We can compare the paths of $\theta$ for the two samplers in side-by-side plots.

//Declare the structure
struct plotControl myMHPlot;

//Fill myMHPlot with 'XY' defaults
myMHPLot = plotGetDefaults("xy");

//Set plot title
plotSetTitle(&myMHPlot, "Metropolis-Hastings Sampler", "Arial", 18);

//Set axis labels
plotSetYLabel(&myMHPlot, "Independent Chain", "Arial", 14);

//Independent chain plot
//Split plot canvas into 2 x 1 grid
//and place next graph in first cell
plotLayout(2, 1, 1);

//Create a sequential vector counting iteration numbers
x_plot = seqa(1, 1, burn_in + keep);

//Plot theta_1
plotXY(myMHPlot, x_plot, theta_draw_ic);

//Random walk plot
//Split plot canvas into 2 x 1 grid
//and place next graph in second cell
plotLayout(2, 1, 2);

//Turn off title for bottom plot
plotSetTitle(&myMHPlot, "");

//Set axis labels
plotSetYLabel(&myMHPlot, "Random Walk");

//Plot theta_1
plotXY(myMHPlot, x_plot, theta_draw_rw);

The above code should produce a graph similar to the image below:

Have a Specific Question?

Get a real answer from a real person

Need Support?

Get help from our friendly experts.

REQUEST A FREE QUOTE

Thank you for your interest in the GAUSS family of products.

© Aptech Systems, Inc. All rights reserved.

Privacy Policy | Sitemap