Introduction to Markov Chain Monte Carlo

Model fitting and inference for infectious disease dynamics

Recap

Recap: Bayesian inference

Last time we saw that the posterior distribution of \(\theta\), given observed data, is

\[p(\theta \mid \text{data}) \propto p(\text{data} \mid \theta)\, p(\theta)\]

Our aim is to draw samples from this distribution.

Why sample?

Sampling to learn about a population

  • Imagine you want to know the mean height of people in London.
  • You cannot measure everyone, so you take a sample \(S\).

\[m = \frac{1}{N_\text{Lon}} \sum_{i \in \text{Lon}} h_i \;\approx\; \frac{1}{N_S} \sum_{i \in S} h_i\]

We can do exactly the same thing to study a posterior distribution.

Monte Carlo sampling

Explore a distribution by drawing from it. More samples, better picture:

Wasting time in sampling

Sampling blindly from the prior spends almost all its effort on parameter values that fit the data badly. We need something that concentrates where the posterior is.

Rejection sampling

The problem

  • Consider a distribution \(f(\theta)\) which we can evaluate for any \(\theta\)
  • How do we draw samples from it?

A proposal distribution

Rejection sampling uses a proposal distribution \(q(\theta)\) which

  • is simple to evaluate
  • is easy to sample from
  • admits some \(M > 1\) with \(f(\theta) < M q(\theta)\) for all \(\theta\)

The algorithm

  1. Sample \(\theta^*\) from \(q(\theta)\)
  2. Draw \(u \sim \text{Uniform}[0,\, M q(\theta^*)]\)
  3. Evaluate \(f(\theta^*)\)
  4. If \(f(\theta^*) > u\) accept, else reject
  5. Repeat steps 1–4

Building up a sample

  1. Sample \(\theta^*\) from \(q(\theta)\)
  2. Draw \(u \sim \text{Uniform}[0,\, M q(\theta^*)]\)
  3. Evaluate \(f(\theta^*)\)
  4. If \(f(\theta^*) > u\) accept, else reject
  5. Repeat steps 1–4

Accepted points build up a sample from \(f(\theta)\).

When rejection sampling struggles

  • Works best when \(q(\theta) \approx f(\theta)\), i.e. \(M \gtrapprox 1\)
  • The acceptance rate is \(1/M\)
  • Requiring \(f(\theta) < M q(\theta)\) everywhere can push the rejection rate very high
  • Worse still in high dimensions

What is Markov chain Monte Carlo?

Markov chain Monte Carlo

  • In Markov chain Monte Carlo (MCMC) we do not need one proposal density \(q(\theta)\) that bounds \(f(\theta)\) everywhere.

  • Instead we build a chain of samples, where each proposed \(\theta^*\) depends on the previous one: the proposal density takes the form \(q(\theta^* \mid \theta)\).

  • A commonly used MCMC algorithm is Metropolis-Hastings.

  • Its acceptance rate is carefully derived to ensure unbiased samples.

Metropolis-Hastings

  1. Initialise \(\theta^{0}\), set \(\theta = \theta^{0}\)
  2. Sample \(\theta^* \sim q(\theta^* \mid \theta)\)
  3. Compute the acceptance probability \(r\)
  4. Draw \(u \sim \text{Uniform}[0,1]\)
  5. Set the new sample to \[\theta^{(s+1)} = \begin{cases} \theta^*, & u < r\\ \theta^{(s)}, & u \geqslant r\end{cases}\]
  6. Repeat steps 2–5

The acceptance probability

If \(q(\theta^* \mid \theta)\) is symmetric, then

\[r = \min\left(1, \frac{f(\theta^*)}{f(\theta)}\right)\]

  • Always move to \(\theta^*\) if it is more probable than \(\theta\)
  • May move even if \(\theta^*\) is less probable

If \(q(\theta^* \mid \theta)\) is asymmetric, then

\[r = \min\left(1, \frac{f(\theta^*)\,q(\theta \mid \theta^*)}{f(\theta)\,q(\theta^* \mid \theta)}\right)\]

Handing it to a library

You will write Metropolis-Hastings once, to know what it does. After that, use a probabilistic programming language.

In Julia that is Turing.jl. You declare the model and it works out the rest:

@model function sir_model(obs)
    R_0 ~ Uniform(1, 20)          # ~ declares a prior
    D_inf ~ Uniform(1, 14)
    traj = simulate(R_0, D_inf)
    obs ~ arraydist(Poisson.(traj))   # ...and a likelihood
end

| conditions the model on data; sample() then runs whichever algorithm you name. The same model definition works with MH, RAM or NUTS.

Better proposals

Tuning by hand does not scale

The proposal covariance has to match the shape of the posterior.

With 2 parameters that is 3 numbers to tune: two variances and a covariance.

With 10 parameters it is 55.

And you generally do not know the posterior scale before you have sampled it.

Adaptive MCMC

  • Adaptive MCMC alters the proposal distribution while the chain is running.

  • Start with a large symmetric variance and scan around to find a mode.

  • Then alter the shape of the proposal to match the covariance of the accepted values.

  • Eventually the proposal density should match the shape of the target density.

Adaptive MCMC: two-stage adaptation

Classical adaptive Metropolis: Haario et al. (2001). Robust adaptive Metropolis (RAM), which we use in the practical: Vihola (2012).

Using gradients

Random walks waste information

Metropolis-Hastings proposes a direction at random and hopes for the best.

But for a model built from operations we can differentiate — including an ODE model, since the solver propagates derivatives — we can compute the gradient of the log-posterior.

The gradient tells us which way is uphill.

Hamiltonian Monte Carlo

Instead of a random step, simulate a puck rolling on the negative log-posterior surface.

  • gradients pull it towards high-density regions
  • momentum carries it a long way before it stops

Both the direction and the length of each move now come from the posterior itself, rather than from a proposal we had to guess.

NUTS

NUTS, the No U-Turn Sampler, is HMC with two things automated:

  • the step size, adapted during warmup
  • the trajectory length, stopped when the path starts turning back on itself — hence the name

The result usually just works, with no proposal tuning at all. It is the default you should reach for.

Hoffman and Gelman (2014); see Betancourt (2017) for the intuition.

Which sampler?

Sampler Best for Limitations
MH Simple models, discrete parameters, learning MCMC Needs manual tuning; slow in many dimensions
RAM Low-dimensional models where gradients are unavailable Still a random walk
NUTS Continuous parameters, differentiable models No discrete parameters; needs gradients

Rule of thumb: start with NUTS. Fall back to RAM or MH when your model has discrete parameters or a stochastic simulator that breaks differentiation — as the particle filter will later in the course.

Your Turn

In the practical you will

  • implement Metropolis-Hastings yourself for a one-dimensional target
  • see what the proposal distribution does to the chain
  • experience how hard tuning a proposal by hand is
  • hand the same job to Turing.jl, then compare MH, RAM and NUTS

References

Betancourt, Michael. 2017. “A Conceptual Introduction to Hamiltonian Monte Carlo.” arXiv Preprint arXiv:1701.02434. https://arxiv.org/abs/1701.02434.
Haario, Heikki, Eero Saksman, and Johanna Tamminen. 2001. “An Adaptive Metropolis Algorithm.” Bernoulli 7 (2): 223–42. https://doi.org/10.2307/3318737.
Hoffman, Matthew D., and Andrew Gelman. 2014. “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.” Journal of Machine Learning Research 15: 1593–623. https://jmlr.org/papers/v15/hoffman14a.html.
Vihola, Matti. 2012. “Robust Adaptive Metropolis Algorithm with Coerced Acceptance Rate.” Statistics and Computing 22 (5): 997–1008. https://doi.org/10.1007/s11222-011-9269-5.

Return to the session