Skip to main content

Markov Chain Monte Carlo (MCMC) refers to a class of simulation methods for generating samples from a complex target distribution by generating random numbers from a Markov Chain whose stationary distribution is the target distribution. MCMC methods are typically used when more direct methods for random number generation (e.g. inversion method) are infeasible. The very first MCMC method was the Metropolis (et al.) algorithm, later expanded by Hastings.

The principle behind Markov Chain Monte Carlo (MCMC) methods is that, by definition, an ergodic recurrent Markov chain $(x_t)_t$ converges in distribution to its stationary distribution no matter what the starting value $x_0$ is. If one can build an ergodic recurrent Markov chain with a given stationary distribution, $\pi$ say, this property can be turned into a simulation principle.

As it happens, there exists a universal algorithm that achieves this goal: it is the Metropolis-Hastings algorithm. Given a target distribution with density $\pi$ (known up to a normalisation constant) and an arbitrary conditional distribution with density $q(y|x)$ on the same space, the Metropolis-Hastings algorithm builds a Markov chain by starting at an arbitrary value $x_0$ (this is the appeal of ergodicity) and by moving from $x_t$ to $x_{t+1}$ according to the rule:

  1. generate a value $y_{t+1}\sim q(y|x_t)$
  2. compute $$ \rho(x_t,y_{t+1})= \text{min} (1, \pi(y_{t+1})q(x_t|y_{t+1})\big/ \pi(x_t)q(y_{t+1}|x_{t}) ) $$
  3. set $$ x_{t+1} = \begin{cases} x_t &\text{with probability }1-\rho(x_t,y_{t+1})\cr y_t &\text{with probability }\rho(x_t,y_{t+1})\cr \end{cases} $$

The validation of this algorithm follows from $\pi$ being stationary and the chain being irreducible if $q$ is everywhere positive.

Another example of MCMC algorithm is the Gibbs sampler. Given a joint distribution $\pi(\theta_1,\ldots,\theta_p)$, moving from $\theta^t$ to $\theta^{t+1}$ proceeds by simulating

  1. $\theta_1^{t+1} \sim \pi(\theta_1|\theta_2^t,\ldots, \theta_p^t)$
  2. ...
  3. $\theta_p^{t+1} \sim \pi(\theta_p|\theta_1^{t+1},\ldots, \theta_{p-1}^{t+1})$

Once again $\pi$ is stationary along this transition.

Two books on the topic by Robert and Casella are Introducing Monte Carlo Methods with R and Monte Carlo Statistical Methods.

More involved versions of MCMC algorithms are the Metropolis-adjusted Langevin algorithm (MALA) and the Hamiltonian Monte Carlo algorithm (HMC), which is behind the Stan software.