Metropolis-Hastings Sampler

GOALS

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

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

INTRODUCTION

The Metroplis-Hastings sampler is an iterative algorithm which produces a sequence of draws θr which can be used to estimate sample parameters such that

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 = 1, . . ., R$ sets θr equal to θ (r-1) , if the draw is rejected.

The Metropolis-Hastings steps

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

Set up sampler specifications, including number of iterations and number of burn-ins draws.
Choose a starting value for θ.
Draw θ $θ^{*}$ from the candidate generating density.
Calculate the acceptance probability α(θ(r-1),θ $α (θ^{(r - 1)}, θ^{*})$ .
Set θ(r) =θ $θ^{(r)} = θ^{*}$ with probability α(θ(r-1),θ $α (θ^{(r - 1)}, θ^{*})$ , or else set θ(r) $θ^{(r)} = θ^{(r - 1)}$ .
Repeat steps 3-5 for r=1 $r = 1, . . ., R$ .
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

The acceptance probability

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

Note : In this example we will set d=6. However, in practice convergence diagnostics should be used to help choose d to optimize the sampler.

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 θ to 1.
Set the standard deviation of the candidate draw to 6.

 
 new;

/*
 ** Set up Metropolis Hastings chain
 ** Discard r0 burnin replications
 */
 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:

Draw a candidate θ using the GAUSS rndn function.
Compute the acceptance probability.
Set θ(r) = θ $θ^{*}$ with probability α(θ(r-1),θ $α (θ^{(r - 1)}, θ^{*})$ , or else set θ(r) = θ(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;

/*
 ** Discard r0 burnin draws
 ** 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

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

The candidate generating density

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

The acceptance probability

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

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;

/*
 ** Discard r0 burnin draws
 ** 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;

Note: It is more efficient to program the random and independent chain sampler in the same for loop. This minimizes the amount of looping in our program. However, we compute a separate for loop for the for the sake of demonstration.

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
 ** Discard r0 burnin replications
 */
 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;

// Discard r0 burn-in draws
 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;

//Discard r0 burn-in draws
 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;

Epidemiologia e Biostatistica

Scienze Sociali

Econometria

ECONOMETRIA FINANZIARIA

Corsi per l'utilizzo del software

Summer school

CONVEGNO ITALIANO DEGLI UTENTI DI STATA

Analisi biostatistica, epidemiologica e ricerca medica

Software per ricerche operative

Analisi statistica generale

formazione multimediale

modelli gerarchici lineari e non lineari

Analisi di data mining

Trasferimento di archivi di dati

Analisi spaziale

Matematica e Ingegneria

word processing scientifico

Analisi statistica specialistica

Disegno di esperimenti e analisi della dimensione dei campioni

Analisi di serie temporali e la stima di modelli econometrici

analisi qualitativa

modelli di reti neurali

STATA PRESS

Altri testi relativi a Stata