The Gibbs sampler requires an ability to generate draws from the conditional distributions. And again, there are many cases where the conditional distributions don’t have convenient Monte Carlo properties. It’s still possible to generate a Markov chain which will converge to the correct distribution, by using techniques similar to importance sampling. The main difference between these is that importance sampling generates independent draws and weights them to achieve the correct distribution. Because the Gibbs sampler is only approximating the distribution at the end of a chain, a weighting system won’t work. Instead, the “weighting” is done by keeping high-density draws for multiple sweeps.

 

While there are quite a few variations of this procedure, there are two principal methods of generating a test draw: you can either draw from a fixed distribution (similar to what is done with importance sampling), or you can draw from a distribution centered around the last draw.

 

The distribution from which you make the draw is known as the proposal distribution or jumping distribution. Let \(x\) be the previous draw and \(y\) be the test draw. You compute a jumping probability \(\alpha\). The “pseudo-code” then would be:

 

if %ranflip(alpha) {

   ...accept y ...

}

else {

   ...keep x ...

}

 

It’s possible for ALPHA to be bigger than 1. If that’s the case, the draw will always be accepted.

 

Fixed proposal distribution (Independence Chain)

Let \(f\) be the target density and \(g\) be the proposal density. (These only need to be known up to a constant of integration). Then

\begin{equation} \alpha = \frac{{f\left( y \right)g\left( x \right)}}{{f\left( x \right)g\left( y \right)}} \end{equation}

As in importance sampling, you need to watch for overflows in the calculation. To safeguard against this, compute this (if possible) as

\begin{equation} \alpha = \exp \left( {\log f\left( y \right) - \log f\left( x \right) + \log g\left( x \right) - \log g\left( y \right)} \right) \end{equation} 

Random Walk Metropolis

If \(y\) is drawn from a density which is conditional on \(x\), let \(g\left( {x,y} \right)\) denote the density for the transition from \(x\) to \(y\). The general formula is

\begin{equation} \alpha = \frac{{f\left( y \right)g\left( {x,y} \right)}}{{f\left( x \right)g\left( {y,x} \right)}} \label{eq:boot_mhgeneralalpha} \end{equation}

However if the transition density has \(y|x\sim N\left( {x,\Sigma } \right)\) where \(\Sigma\) doesn’t depend upon \(x\), then \(g\left( {y,x} \right) = g\left( {x,y} \right)\), so \eqref{eq:boot_mhgeneralalpha} simplifies to

\begin{equation} \alpha = \frac{{f\left( y \right)}}{{f\left( x \right)}} \end{equation}

This is known as Random Walk Metropolis. This simplification will happen whenever \(g\left( {y,x} \right) = g\left( {x,y} \right)\), which generally would mean a Normal or t centered at the previous draw. Note that if you do any linearization to get the proposal density, that will make the density dependent upon \(x\) for more than the mean. If you do that, you have to do the analogous linearization around \(y\) to be able to compute \(g\left( {y,x} \right)\).

 

Choosing a Proposal Density

This can be quite difficult, and may require quite a bit of experimentation, particularly if you’re applying this to a multi-dimensional \(x\). With either type of proposal, if the distribution is too broad, you don’t move very often because most of the draws hit low densities. If it’s too narrow, you might move only very slowly once you’ve landed in a high density area. It’s a good idea to keep track of the percentage of times you jump. If this is too small, it may indicate that you need to make some type of adjustment. But taking many jumps doesn’t necessarily indicate a successful set-up—you can move around quite a bit in one area, but fail to move beyond it.

 

Displaying the Acceptance Rate

If you are doing a time-consuming chain, it’s helpful to get quicker feedback on the chain’s acceptance probability than you would get by running an entire set of draws first. You can do that by replacing the INFOBOX inside the loop with:

 

infobox(current=draw) %strval(100.0*accept/(draw+nburn+1),"##.#")

 

This puts a second line in the progress box showing the percentage of jumps made. This is updated after each sweep, so if it looks like it’s running too low, you can cancel the loop, and try some different settings.

Examples

For a (relatively simple) example, see GARCHGIBBS.RPF. MONTESVAR_MH.RPF analyzing a structural VAR using Metropolis-Hastings.

Copyright © 2025 Thomas A. Doan