### Introduction

In this tutorial, we will

- Implement the Metropolis-Hastings sampler using two different variations of acceptance probabilities:
- An independence chain
- A random walk chain

- Compare the samples from the two methods in side-by-side histograms.
- 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:

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

- Draw a candidate $\theta$ using the GAUSS
`rndn`

function. - Compute the acceptance probability.
- 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;
```

## Print results

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: