Markov chain Monte Carlo (MCMC) sampling

Estimated time: 115 minutes

Slides

Introduction

In the previous session, you learned Bayesian inference through manual calculations: computing priors, likelihoods, and posteriors and exploring the posterior distribution using grid search. This approach introduced you to the core components of conducting inference, but it also exposed a fundamental problem: grid search doesn’t scale.

This session introduces Markov Chain Monte Carlo (MCMC), the algorithmic solution that makes Bayesian inference practical for real problems. Rather than exhaustively searching parameter space, MCMC intelligently explores it by taking random walks that spend more time in high-probability regions.

You’ll build understanding by walking through a simple MCMC implementation, then see how modern tools like Turing.jl automate everything with state-of-the-art methods. By the end, you’ll understand both what MCMC does under the hood and how to use it effectively in practice.

Objectives

The aim of this session is to learn how to sample from posterior distributions using MCMC. By building the algorithm yourself first, you’ll understand what modern tools like Turing.jl do under the hood.

In this session, you will:

  1. Implement the Metropolis-Hastings MCMC algorithm and apply it to a simple problem
  2. Experience the challenges of using MCMC: proposal tuning and mixing
  3. See how Turing.jl automates Metropolis-Hastings for you
  4. Learn why NUTS is the modern upgrade that solves some of the problems with Metropolis-Hastings
  5. Apply NUTS to conduct Bayesian inference on the SIR model

Setup

Source file

The source file of this session is located at sessions/mcmc.qmd.

Load libraries

using Random ## for random numbers
using Distributions ## for probability distributions
using DifferentialEquations ## for differential equations
using DataFrames ## for data frames
using Plots ## for plots
using StatsPlots ## for statistical plots
using Turing ## for probabilistic programming and MCMC
using FlexiChains ## chain objects and diagnostics
using AdvancedMH ## for adaptive MCMC samplers

Initialisation

We set a random seed for reproducibility. Setting this ensures that you should get exactly the same results on your computer as we do.

Random.seed!(1234)
TaskLocalRNG()

Load data and model from previous session

We’ll use the same SIR model and epidemic data from the introduction session.

Key functions:

  • sir_ode!(du, u, p, t) - ODE system for S, I, R compartments
  • simulate_sir(θ, init_state, times) - Returns DataFrame with trajectory
using DifferentialEquations
using Distributions
using DataFrames

"""
    sir_ode!(du, u, p, t)

SIR model ordinary differential equations.

# States
- `u[1]` (S): Susceptible
- `u[2]` (I): Infectious
- `u[3]` (R): Recovered

# Parameters
- `p[1]` (R_0): Basic reproduction number
- `p[2]` (D_inf): Infectious period (days)
"""
function sir_ode!(du, u, p, t)
    S, I, R = u
    R_0, D_inf = p
    N = S + I + R

    β = R_0 / D_inf
    ν = 1 / D_inf

    du[1] = -β * S * I / N          # dS/dt
    du[2] = β * S * I / N - ν * I   # dI/dt
    du[3] = ν * I                   # dR/dt
end

"""
    SIR_BASE_PROBLEM

A template `ODEProblem` for `sir_ode!`, built once when the package loads.
`simulate_sir` calls `remake` on it rather than building a fresh problem on
every call, which matters because a likelihood asks for tens of thousands of
simulations. The placeholder state and parameters are never used: `remake`
replaces `u0`, `p` and `tspan` on every call.
"""
const SIR_BASE_PROBLEM = ODEProblem(sir_ode!, zeros(3), (0.0, 1.0), zeros(2))

"""
    simulate_sir(θ, init_state, times)

Simulate the deterministic SIR model.

# Arguments
- `θ`: Dict with keys :R_0, :D_inf
- `init_state`: Dict with keys :S, :I, :R
- `times`: Time points to return (e.g., 0.0:1.0:30.0)

# Returns
DataFrame with columns: time, S, I, R, Inc (daily incidence)

`Inc` is the flow out of S, so `Inc[i]` counts the infections between
`times[i - 1]` and `times[i]`. No such interval precedes the first time point,
so `Inc[1]` is zero and `times` needs a point before the first observation:
`Inc[i + 1]` is the incidence over the day the i-th observation counts.
"""
function simulate_sir(θ, init_state, times)
    u0 = [init_state[:S], init_state[:I], init_state[:R]]
    prob = remake(
        SIR_BASE_PROBLEM;
        u0 = u0,
        p = [θ[:R_0], θ[:D_inf]],
        tspan = (times[1], times[end]),
    )
    sol = solve(prob, Tsit5(), saveat = times)

    df = DataFrame(
        time = sol.t,
        S = [u[1] for u in sol.u],
        I = [u[2] for u in sol.u],
        R = [u[3] for u in sol.u],
    )

    # Daily incidence: new infections, the flow out of S
    df.Inc = [0.0; -diff(df.S)]

    return df
end
# Initial state and time range
init_state = [999, 1, 0]
t0, tend = (0, 30)

# Epidemic data
epi1 = DataFrame(
    time = 1:1:30,
    obs = [
        1, 3, 4, 11, 19, 38, 67, 112, 183, 210, 237, 182, 178, 154, 95, 68, 51,
        30, 32, 18, 19, 9, 2, 2, 2, 0, 0, 1, 0, 0
    ]
)
30×2 DataFrame
5 rows omitted
Row time obs
Int64 Int64
1 1 1
2 2 3
3 3 4
4 4 11
5 5 19
6 6 38
7 7 67
8 8 112
9 9 183
10 10 210
11 11 237
12 12 182
13 13 178
⋮ ⋮ ⋮
19 19 32
20 20 18
21 21 19
22 22 9
23 23 2
24 24 2
25 25 2
26 26 0
27 27 0
28 28 1
29 29 0
30 30 0

From manual exploration to algorithmic sampling

In the previous session, we manually explored the posterior by trying different parameter combinations. This worked for our simple SIR model with 2 parameters, but we also saw that the number of combinations to explore can be prohibitive.

Instead of trying all combinations, Markov chain Monte Carlo works in the following way:

  1. Start somewhere in parameter space
  2. Propose a small move to nearby parameters
  3. Accept or reject the move based on how much better or worse the posterior is
  4. Repeat many times to build up a sample

This is the Metropolis-Hastings algorithm.

Metropolis-Hastings for a simple case

At its heart, Metropolis-Hastings is a method for drawing samples from any probability distribution we can evaluate at every point, even when we can’t sample from it directly. That includes posteriors, but also familiar distributions like the Gaussian. Before we point it at anything complicated, let’s run it on Normal(5, 2) so we can compare the samples against a known truth:

# Target distribution: Normal(5, 2)
target_mean = 5.0
target_std = 2.0
target_dist = Normal(target_mean, target_std)

# Helper function to evaluate log-density of target
function eval_log_density(x)
    return logpdf(target_dist, x)
end

# Simple Metropolis-Hastings sampler
function metropolis_hastings_1d(log_density_fn, n_samples, initial_value,
                                proposal_std)
    # Initialise sample vector
    samples = Float64[]
    # current_x tracks the state of the chain
    # initially, set to given initial value
    current_x = initial_value
    # Evaluate log probability at the initial value
    current_log_prob = log_density_fn(current_x)
    # Track number of acceptances
    n_accepted = 0

    for i in 1:n_samples
        # Propose new state
        proposal_x = current_x + randn() * proposal_std
        # Evaluate log probability at target
        proposal_log_prob = log_density_fn(proposal_x)

        # Calculate acceptance probability (in log space)
        # Note: This is the Metropolis ratio (symmetric proposal)
        # Full Metropolis-Hastings would include proposal ratio q(x|x')/q(x'|x)
        # But Normal proposal is symmetric, so this ratio = 1 and cancels out
        log_alpha = proposal_log_prob - current_log_prob
        alpha = exp(log_alpha)

        # Accept or reject
        if rand() < alpha
            current_x = proposal_x
            current_log_prob = proposal_log_prob
            n_accepted += 1
        end # if reject, current_x stays what it is

        push!(samples, current_x)
    end

    acceptance_rate = n_accepted / n_samples
    return samples, acceptance_rate
end

# Test our sampler
samples_mh, acceptance_rate = metropolis_hastings_1d(
    eval_log_density, 5000, 0.0, 1.0
)

println("Acceptance rate: $(round(acceptance_rate * 100, digits = 1))%")
println("Sample mean: $(round(mean(samples_mh), digits = 2)) (true: $(target_mean))")
println("Sample std: $(round(std(samples_mh), digits = 2)) (true: $(target_std))")
Acceptance rate: 84.0%
Sample mean: 5.12 (true: 5.0)
Sample std: 2.08 (true: 2.0)
# Visualise the results
p1 = plot(
    samples_mh, title = "Trace plot", xlabel = "Iteration", ylabel = "x",
    linewidth = 1, label = ""
)
hline!(p1, [target_mean], color = :red, linewidth = 2, label = "True mean")

p2 = histogram(
    samples_mh, normalize = :pdf, alpha = 0.7,
    title = "Sampled distribution", xlabel = "x", ylabel = "Density",
    label = "MCMC samples"
)

# Overlay true distribution
x_range = range(-5, 15, length = 100)
plot!(
    p2, x_range, pdf.(target_dist, x_range), linewidth = 3, color = :red,
    label = "True distribution"
)

plot(p1, p2, layout = (2, 1), size = (800, 600))

The trace plot shows the chain rapidly settling into the region around the true mean and then bouncing around it — a well-mixed chain is often described as looking like a fat, hairy caterpillar that fills the vertical band around its target. The histogram of samples closely matches the true density.

Notice the acceptance rate, which is over 80%. That is comfortable rather than efficient: the chain accepts nearly everything it proposes because it proposes such small steps. Widening the proposal buys bigger steps at the price of more rejections, and the section on proposal distributions below works out where the best trade lies.

TipExercise: What is the proposal distribution?

Look back at metropolis_hastings_1d. Which line draws the proposed new state, and what distribution does it use?

The line proposal_x = current_x + randn() * proposal_std draws a proposed new state. randn() returns a standard Normal sample, so current_x + randn() * proposal_std is a draw from Normal(current_x, proposal_std) — the proposal distribution is a Normal centred on the current state. proposal_std controls how far the sampler tries to jump on each step. The next subsection explores what happens when we make it too small or too large.

The Julia MH sampler above translates directly:

target_mean <- 5.0
target_std <- 2.0

eval_log_density <- function(x) {
  dnorm(x, mean = target_mean, sd = target_std, log = TRUE)
}

metropolis_hastings_1d <- function(log_density_fn, n_samples,
                                   initial_value, proposal_std) {
  samples <- numeric(n_samples)
  current_x <- initial_value
  current_log_prob <- log_density_fn(current_x)
  n_accepted <- 0

  for (i in 1:n_samples) {
    proposal_x <- current_x + rnorm(1, mean = 0, sd = proposal_std)
    proposal_log_prob <- log_density_fn(proposal_x)

    log_alpha <- proposal_log_prob - current_log_prob
    if (runif(1) < exp(log_alpha)) {
      current_x <- proposal_x
      current_log_prob <- proposal_log_prob
      n_accepted <- n_accepted + 1
    }
    samples[i] <- current_x
  }

  list(samples = samples, acceptance_rate = n_accepted / n_samples)
}

result <- metropolis_hastings_1d(eval_log_density, 5000, 0.0, 1.0)

The general principle: MCMC as a universal random number generator

What we just built solves a fundamental problem in computational statistics: given a probability distribution where we can evaluate the density at any point, how do we generate random samples from it?

Many useful probability distributions don’t have simple sampling algorithms. Easy distributions like Normal (randn()) and Uniform (rand()) have built-in samplers. But posteriors, mixtures, and custom models have no simple formula.

MCMC solves this by only requiring that we can evaluate the log-probability at any point. We used a Normal distribution above (even though we could just use randn()!) to demonstrate this general principle.

Even “simple” distributions like the Normal don’t have a direct sampling formula. When you call randn() in Julia or rnorm() in R, the computer uses transformation methods:

Common approaches:

  1. Box-Muller transform: Convert two uniform random numbers to two normal random numbers

    # Simplified version
u1, u2 = rand(), rand()
z1 = sqrt(-2 * log(u1)) * cos(2π * u2)
z2 = sqrt(-2 * log(u1)) * sin(2π * u2)
# z1 and z2 are standard normal

  2. Inverse CDF method: Transform uniform → normal via inverse cumulative distribution function

  3. Ziggurat algorithm: Fast, more complex method (what Julia actually uses)

Even “easy” distributions require clever tricks. For complex posteriors with correlations, ODEs, and constraints, no such formula exists, and we need a method like MCMC.

The posterior distribution is just another probability distribution. We can evaluate it at any point: log_posterior(θ) = log_prior(θ) + log_likelihood(θ|data). Because of this, MCMC can sample from it. The samples give us everything we need: visualisation (plot the distribution shape), summaries (mean, median, credible intervals), and forward simulation (with parameter correlations preserved).

When we apply MCMC to the SIR model later, we will be doing exactly the same thing as we did here for the Normal example - just with a more complex probability distribution that happens to be a Bayesian posterior. MCMC doesn’t “know” about Bayesian inference; it just knows how to sample from distributions where we can evaluate the density.

Understanding proposal distributions

The choice of proposal standard deviation matters:

# Test different proposal standard deviations
proposal_stds = [0.1, 0.5, 1.0, 5.0, 10.0, 50.0]
results = []

for prop_std in proposal_stds
    samples, acc_rate = metropolis_hastings_1d(eval_log_density, 5000, 0.0, prop_std)
    push!(results, (std = prop_std, samples = samples, acc_rate = acc_rate))
end

# Print acceptance rates
println("Acceptance rates for different proposal std:")
for result in results
    println(
        "  Proposal std = $(result.std): ",
        "$(round(result.acc_rate*100, digits = 1))%"
    )
end

# Plot trace plots for different proposal stds
plots = Plots.Plot[]
for (i, result) in enumerate(results)
    p = plot(
        result.samples,
        title = "Prop. std = $(result.std), " *
            "Acc = $(round(result.acc_rate*100, digits = 1))%",
        xlabel = "Iteration", ylabel = "x", linewidth = 1, ylim = (-5, 15)
    )
    hline!(p, [target_mean], color = :red, linewidth = 2, label = "")
    push!(plots, p)
end

plot(plots..., layout = (length(plots), 1), size = (800, 1000))
Acceptance rates for different proposal std:
  Proposal std = 0.1: 98.5%
  Proposal std = 0.5: 91.8%
  Proposal std = 1.0: 85.1%
  Proposal std = 5.0: 42.6%
  Proposal std = 10.0: 24.2%
  Proposal std = 50.0: 5.0%
ImportantChoosing the proposal width
  • Too small proposals give high acceptance but slow exploration
  • Too large proposals give low acceptance and wasted computation
  • The best acceptance rate depends on how many parameters you are sampling: around 44% for a single parameter, falling towards 23% as the number grows

The runs above show this. The chain mixes best at a proposal standard deviation of 5, which accepts around 43% of what it proposes, while the run at a standard deviation of 10 accepts 24% and wanders visibly more slowly. The 23% figure is the limit for a large number of parameters, so it is the wrong target for the single parameter sampled here.

For the normal distribution, sampling was straightforward. But efficient sampling from complex posterior distributions requires a good proposal distribution, which can be difficult to find. This is where more sophisticated algorithms that tune the proposal distribution automatically to a given problem really help.

NoteGood time for a break

You’ve now seen Metropolis-Hastings from scratch and understand its core mechanics. The next sections show how modern tools automate this. If you’re working through this in one sitting, this is a natural pause point.

Automating MCMC with Turing.jl

Now that you understand how MCMC works internally, let’s see how Turing.jl automates everything you just coded.

Turing.jl on the Normal distribution

First, let’s use Turing.jl to sample from the exact same Normal distribution you just coded manually:

# Define the same Normal(5, 2) distribution in Turing
@model function normal_model()
    x ~ Normal(5.0, 2.0)
end
  • @model - a Turing.jl macro that defines a probabilistic model. Everything inside becomes part of the model specification.
  • x ~ Normal(5.0, 2.0) - the ~ operator means “is distributed as”. This declares that x is a random variable drawn from a Normal distribution with mean 5 and standard deviation 2. Turing will infer the value of x during sampling.

If you’ve used probabilistic programming languages, you’ll recognise the ~ syntax. Turing’s @model is similar to model definitions in Stan, NIMBLE, and BUGS/JAGS:

# NIMBLE / BUGS / JAGS
model({
  x ~ dnorm(5, sd = 2)
})
// Stan
model {
  x ~ normal(5, 2);
}

The key advantage of Turing is that the model is written in native Julia, so you can use loops, conditionals, and call any Julia function (like our ODE solver) directly inside the model — no separate modelling language needed.

Now let’s run the same Metropolis-Hastings algorithm as before, but with Turing handling the plumbing:

# Sample using Metropolis-Hastings, with the same proposal width and the same
# number of iterations as the sampler we wrote by hand
chain_normal_mh = sample(normal_model(), MH([1.0;;]), 5000)
  • sample(model, sampler, n) - runs the MCMC sampler on the model for n iterations, returning a chain of samples.
  • MH(proposal) - Metropolis-Hastings sampler. MH takes the proposal covariance. The [1.0;;] syntax creates a 1×1 matrix, where the double semicolon ;; tells Julia to build a matrix rather than a vector. The value sets the proposal variance, and since the variance is 1.0 here it happens to equal the standard deviation of 1.0 we passed to metropolis_hastings_1d. That coincidence is only true at 1.0: a variance of 0.04 is a standard deviation of 0.2.
# Compare to your manual implementation
println("Manual M-H:")
println("  Mean: $(round(mean(samples_mh), digits = 2))")
println("  Std: $(round(std(samples_mh), digits = 2))")
println("\nTuring.jl M-H:")
summarystats(chain_normal_mh)
Manual M-H:
  Mean: 5.12
  Std: 2.08

Turing.jl M-H:
╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (1) ── AbstractPPL.VarName                                        │
│  Float64  x                                                                  │
│                                                                              │
│ Extras (4)                                                                   │
│  Float64  accepted, logprior, loglikelihood, logjoint                        │
│                                                                              │
│ Summary                                                                      │
│   param    mean     std    mcse  ess_bulk  ess_tail    rhat      q5  …       │
│       x  5.1728  2.0112  0.1588  161.5821  435.2782  1.0077  1.7680  …       │
╰──────────────────────────────────────────────────────────────────────────────╯
  • summarystats(chain) - returns a table with the mean, the standard deviation, MCSE (Monte Carlo standard error), ESS (effective sample size), R̂ (a convergence diagnostic) and the 5%, 50% and 95% quantiles, for every parameter in the chain. The MCMC diagnostics session explains what each of these means and when to worry about it.
  • chain[:x] - extracts the samples for parameter x from the chain. You can use this to access any parameter defined in your @model function.

If you’ve used Stan, the workflow is similar: sample() runs inference and returns a chain object, summarystats() gives summary statistics, and chain[:x] extracts samples for a parameter.

Action rstan cmdstanr Turing.jl
Run sampler sampling(model, ...) model$sample(...) sample(model, NUTS(), n)
Summary summary(fit) fit$summary() summarystats(chain)
Extract parameter as.array(fit, pars="x") fit$draws("x") chain[:x]
Plotting code 
# Visualise - same result!
p1 = histogram(
    samples_mh, normalize = :pdf, alpha = 0.7,
    title = "Manual M-H", xlabel = "x", ylabel = "Density", label = "Manual"
)
p2 = histogram(
    chain_normal_mh[:x], normalize = :pdf, alpha = 0.7,
    title = "Turing M-H", xlabel = "x", ylabel = "Density", label = "Turing"
)
plot(p1, p2, layout = (1, 2), size = (800, 300))

Turing.jl ran the same algorithm with the same proposal width, in a single call to the sample() function. The two chains start from different points, because Turing draws its starting value from the model rather than taking the 0.0 we passed by hand, so the samples are not identical draw for draw. The distributions they explore are the same.

Applying Turing.jl to the SIR model

Turing.jl uses the @model macro to define Bayesian models. This is the same structure we’ve been using manually (priors, simulation, likelihood) but now in a form that Turing can automatically sample from:

sir_prob = ODEProblem(sir_ode!, [999, 1, 0], (0.0, 1.0), [2.5, 5.0])

# Define our SIR model in Turing.jl
@model function sir_model(init_state, times, n_obs)
    # Weakly informative priors
    R_0 ~ truncated(Normal(2.5, 1.5), lower=1.0)
    D_inf ~ truncated(Normal(5.0, 3.0), lower=0.5)

    prob = remake(
        sir_prob, u0 = init_state, p = [R_0, D_inf], tspan = (times[1], times[end])
    )
    sol = solve(prob, saveat = times)

    # Likelihood (Poisson observation model of prevalence)
    I = sol[2, :]
    # max(..., 0.0) guards against tiny negative values from the ODE solver
    # (numerical noise can push I slightly below zero near the end of an epidemic)
    lambdas = [max(I[i], 0.0) for i in 1:n_obs]
    obs ~ arraydist(Poisson.(lambdas))
end

# Create model instance
model = sir_model(init_state, epi1.time, length(epi1.obs)) | (; obs = epi1.obs);

ODEProblem does more than store the arguments you hand it. It wraps the right hand side in an ODEFunction, settles the problem type and builds the indexing machinery the solver interface expects. For a system this small that setup costs about as much as the solve.

remake takes an existing problem and returns a copy with some fields replaced. It does not modify the one you pass it, so a single base problem is safe to share across threads. Declaring it const makes no measurable difference: it is read once per call, immediately before a solve that costs several microseconds.

Two new pieces of syntax appear here:

obs ~ arraydist(Poisson.(lambdas)) — declares a vector of observations where each element has its own distribution. arraydist creates an independent product distribution: obs[1] ~ Poisson(lambdas[1]), obs[2] ~ Poisson(lambdas[2]), etc. The dot in Poisson.(lambdas) is Julia’s broadcasting — it applies Poisson to each element of lambdas.

sir_model(...) | (; obs = epi1.obs) — the | operator conditions the model on observed data. It tells Turing “don’t sample obs — use these observed values instead and compute the likelihood”. The (; obs = epi1.obs) is a named tuple matching the variable name obs inside the model to the data — the semicolon is Julia’s syntax for a named tuple with no positional arguments.

This separation of model structure from data is a key Turing pattern: the @model function defines the generative process, and | supplies the observations. The same model can be used unconditioned (for simulation), conditioned (for inference), or deconditioned after fitting for posterior predictive checks (decondition() — you’ll see this in the Model Checking session).

Sample the SIR model with Metropolis-Hastings

Now let’s use the Metropolis-Hastings algorithm on our epidemiological model:

# Same Metropolis-Hastings algorithm, but automated
chain_mh = sample(model, MH([1.0 0.0; 0.0 1.0]), 5000)
  • [1.0 0.0; 0.0 1.0] - a 2×2 matrix in Julia. Spaces separate columns, semicolons separate rows. This is the identity matrix, used here as the proposal covariance for MH. The diagonal values are variances, one per parameter, so 1.0 here is a proposal standard deviation of 1.0; off-diagonal zeros mean proposals are independent.
println("Turing.jl Metropolis-Hastings results:")
summarystats(chain_mh)
Turing.jl Metropolis-Hastings results:
╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (2) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_inf                                                         │
│                                                                              │
│ Extras (4)                                                                   │
│  Float64  accepted, logprior, loglikelihood, logjoint                        │
│                                                                              │
│ Summary                                                                      │
│   param    mean     std    mcse  ess_bulk  ess_tail    rhat      q5  …       │
│     R_0  2.4263  0.0957  0.0068   18.3776   12.1986  1.2960  2.4065  …       │
│   D_inf  1.9286  0.1729  0.0059   17.9370   12.2195  1.1903  1.9095  …       │
╰──────────────────────────────────────────────────────────────────────────────╯
# Visualise
plot(
    chain_mh[:R_0], title = "Turing MH: R_0 Trace", xlabel = "Iteration",
    ylabel = "R_0", linewidth = 1
)

The trace looks poor: the chain gets stuck at single values for long stretches instead of exploring freely. Compare this against the healthy trace from the Normal example above; that’s the target shape. The identity proposal we passed is a bad match for this posterior, and we can do better by tuning the proposal covariance.

TipExercise: Tune the proposal distribution

Your task is to improve the mixing above by adjusting the proposal covariance matrix. Finding good proposals by hand is hard; the goal is to experience the challenge, not to find a perfect solution.

Suggested approach:

  • Spend 5–10 minutes experimenting with different proposal matrices.
  • If you’re stuck after 10 minutes, reveal the solution — you’ve learned the key lesson.
  • The frustration is the point: manual tuning doesn’t scale to real problems.

The proposal covariance matrix controls how the sampler proposes new parameter values:

# Current: identity matrix (proposal std = 1.0 for both parameters)
MH([1.0 0.0; 0.0 1.0])

# Structure: [var_R0  cov; cov  var_Dinf]
# Diagonal = variance (std²) for each parameter
# Off-diagonal = covariance between parameters

Try these experiments (run each for 5000 samples):

  1. Smaller proposals: MH([0.01 0.0; 0.0 0.01]) — what happens to acceptance rate?
  2. Larger proposals: MH([5.0 0.0; 0.0 5.0]) — does it explore faster or slower?
  3. Even smaller: MH([0.001 0.0; 0.0 0.001]) — the posterior turns out to be very concentrated for this dataset.
  4. Add correlation: MH([0.001 0.0009; 0.0009 0.001]) — if parameters are correlated in the posterior, correlated proposals may help.

Questions to answer:

  • What acceptance rate gives the best mixing? (Hint: with a proposal shaped like the posterior, two parameters do best near 35%; with the diagonal proposals above, nearer 20–25%.)
  • Can you find settings that give a well-mixed trace for both parameters simultaneously?
  • How long did tuning take you? Imagine doing this for 10+ parameters.

Acceptance rate too high (>50%)? Your proposals are too small. The chain accepts most moves but takes tiny steps; you’ll see the trace plot moving very slowly.

Acceptance rate too low (<10%)? Your proposals are too large. Most proposed jumps land in low-probability regions and get rejected; you’ll see the trace plot stuck for many iterations.

Sweet spot: a healthy balance of exploration and acceptance. Where it falls depends on the shape of the proposal as well as the number of parameters: around 20–25% for the diagonal proposals above, and nearer 35–45% once the proposal is aligned with the posterior.

Check your acceptance rate with: mean(diff(vec(chain[:R_0])) .!= 0)

The identity proposal we started with uses std = 1.0 for each parameter, but this posterior is very concentrated (both parameters have std around 0.04 once you get there). Proposal std much larger than the posterior width means almost every jump lands in a low-density region and gets rejected.

Rule of thumb: proposal std ≈ 0.5–1× posterior std, i.e. proposal variance around 0.001 for each parameter here.

Run the chain and look at a scatter plot of R₀ vs D_inf samples. Here they are strongly positively correlated (higher R₀ → higher D_inf, correlation about 0.9), so independent proposals are inefficient.

Imagine the posterior is a narrow diagonal ridge. Independent proposals move horizontally or vertically, often stepping off the ridge. Correlated proposals can move along the ridge.

Try an off-diagonal value of the same sign as the correlation, e.g. MH([0.001 0.0009; 0.0009 0.001]).

# Space for your experiments - modify the proposal matrix and re-run
# chain_tuned = sample(model, MH([??? ???; ??? ???]), 5000)
# plot(chain_tuned)

After experimentation, a proposal like this works reasonably well:

# Proposal variance ~0.001, with a strong positive correlation
chain_tuned = sample(model, MH([0.001 0.0009; 0.0009 0.001]), 5000)

Why this works:

  1. Small proposals (variance 0.001, std ≈ 0.03): the posterior for both parameters has standard deviation ≈ 0.04 for this dataset — much narrower than the priors suggest. Any proposal much wider than that gets rejected almost every time.
  2. Positive correlation: the data pin down the epidemic’s early growth rate, roughly (R0−1)/Dinf(R_0 - 1) / D_{inf}, so a longer infectious period needs a larger R₀ to produce the same epidemic. The posterior is therefore a ridge running along the positive diagonal, with a correlation of about 0.9, and a proposal aligned with it moves along the ridge instead of stepping off it. On this dataset that is worth about seven times the effective sample size of the negatively correlated proposal, and about three times that of the plain diagonal one.

The difficulty you just experienced is why manual tuning doesn’t scale. With 2 parameters, you have 3 values to tune (two variances, one covariance). With 10 parameters, you’d have 55 values. This motivates adaptive methods that learn the proposal automatically.

Adaptive MCMC: learning better proposals

The problem with standard MH is that we need to manually tune the proposal covariance. If parameters are correlated, independent proposals (diagonal covariance) will mix poorly.

Adaptive Metropolis methods sidestep this by learning the proposal covariance from the samples themselves as the chain runs. The classical version (Haario et al., 2001) uses the running empirical covariance of the samples seen so far. Robust Adaptive Metropolis (RAM; Vihola, 2012) instead scales the proposal to hit a target acceptance rate, which handles heavy-tailed posteriors better.

Let’s run it on our SIR model:

using LinearAlgebra
# Adaptive MCMC - learns proposal covariance automatically.
# Seed the adaptation with a reasonable initial scale so it doesn't spend
# thousands of steps trying to find the posterior region, and give it 5000
# warmup iterations to adapt in. Those are discarded, so the returned chain is
# the 5000 iterations that follow, sampled with the adapted proposal.
initial_scale = 0.05 .* Matrix{Float64}(I, 2, 2)
chain_ram = sample(
    model,
    externalsampler(AdvancedMH.RobustAdaptiveMetropolis(S = initial_scale)),
    5000;
    num_warmup = 5000,
    check_model = false
)
  • RobustAdaptiveMetropolis(S = ...) - from AdvancedMH.jl. S is a lower-triangular Cholesky factor of the initial proposal covariance; the sampler updates it during warmup, growing or shrinking it to hit a target acceptance rate, and then holds it fixed for the sampling that follows.
  • externalsampler() - wraps an AdvancedMH.jl sampler for use with Turing models, passing the warmup phase through so that RAM has iterations to adapt in.
  • num_warmup = 5000 - 5000 adaptation iterations. They are discarded, so the returned chain is the 5000 that follow, drawn with the adapted proposal.
  • check_model = false - skips model validation (needed for some external samplers).
# Compare MH vs Adaptive MH
p1 = plot(chain_mh[:R_0], title="Standard MH: R_0", xlabel="Iteration", ylabel="R_0", legend=false)
p2 = plot(chain_ram[:R_0], title="Adaptive MH (RAM): R_0", xlabel="Iteration", ylabel="R_0", legend=false)
plot(p1, p2, layout=(1, 2), size=(800, 300))

Notice how the adaptive sampler explores the parameter space more efficiently.

NUTS: a gradient-based upgrade

Even with adaptation, random-walk proposals are limited: they don’t use any information about which direction the posterior is denser in. For any model built from operations we can differentiate (which includes our ODE model, since ODE solvers propagate derivatives) we can compute the gradient of the log-posterior with respect to the parameters. The gradient points from the current parameter values toward regions of higher posterior density.

Hamiltonian Monte Carlo (HMC) uses those gradients to construct proposals that move along the shape of the posterior instead of jumping around blindly. Intuitively, it simulates a puck rolling on the negative log-posterior surface: gradients pull it into high-density regions and momentum lets it travel a long way before proposing a new state. This solves the two big problems of random-walk MH (poorly aligned proposals and short steps) because both the direction and the length of each move come from the posterior itself.

NUTS (the No U-Turn Sampler) is the algorithm we actually run. It’s HMC with two things automated: the step size (adapted during warmup) and the trajectory length (stopped when the simulated puck starts turning back on itself, hence the name). The result is a sampler that usually just works, with no tuning to speak of:

# Modern gradient-based MCMC
# Expect: ~30-60 seconds on first run (includes compilation), ~10-20 seconds on subsequent runs
chain_nuts = sample(model, NUTS(), 2000)
  • NUTS() - the No U-Turn Sampler, a gradient-based MCMC algorithm that automatically tunes its step size and trajectory length. Much more efficient than Metropolis-Hastings for continuous parameters. The algorithm uses automatic differentiation to compute gradients of the log-posterior.

How do gradients help with sampling? Metropolis-Hastings proposes random steps and hopes for the best — like exploring a dark room by stumbling around. NUTS computes the gradient (slope) of the log-posterior, which tells it which direction leads to higher probability. It then simulates a puck sliding along this surface, following the natural contours of the posterior. This lets it take large, efficient steps along ridges and valleys instead of the small, random steps of MH. The “No U-Turn” part automatically decides when to stop sliding, so you don’t need to tune a trajectory length.

println("NUTS results:")
summarystats(chain_nuts)
NUTS results:
╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (2) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_inf                                                         │
│                                                                              │
│ Extras (14)                                                                  │
│  Float64  n_steps, is_accept, acceptance_rate, log_density,                  │
│           hamiltonian_energy, hamiltonian_energy_error,                      │
│           max_hamiltonian_energy_error, tree_depth, numerical_error,         │
│           step_size, nom_step_size, logprior, loglikelihood, logjoint        │
│                                                                              │
│ Summary                                                                      │
│   param    mean     std    mcse  ess_bulk  ess_tail    rhat      q5  …       │
│     R_0  2.4456  0.0365  0.0016  540.6291  524.9272  1.0132  2.3882  …       │
│   D_inf  1.9446  0.0360  0.0015  551.5183  516.3810  1.0117  1.8876  …       │
╰──────────────────────────────────────────────────────────────────────────────╯
Plotting code 
# Compare Metropolis-Hastings vs NUTS
p1 = plot(
    chain_mh[:R_0], title = "Metropolis-Hastings: R_0", xlabel = "Iteration",
    ylabel = "R_0", linewidth = 1, legend = false
)
p2 = plot(
    chain_nuts[:R_0], title = "NUTS: R_0", xlabel = "Iteration",
    ylabel = "R_0", linewidth = 1, legend = false
)

p3 = plot(
    chain_mh[:D_inf], title = "Metropolis-Hastings: D_inf",
    xlabel = "Iteration", ylabel = "D_inf", linewidth = 1, legend = false
)
p4 = plot(
    chain_nuts[:D_inf], title = "NUTS: D_inf", xlabel = "Iteration",
    ylabel = "D_inf", linewidth = 1, legend = false
)

plot(p1, p2, p3, p4, layout = (2, 2), size = (800, 600))

Notice that:

  • NUTS mixes better: it explores the parameter space more efficiently.
  • NUTS needs fewer samples for the same quality.
  • NUTS had automatic warmup: no manual burn-in needed, see below.

MCMC chains start at an arbitrary point in parameter space, which may be far from the high-probability region. The initial samples as the chain “finds its way” to the posterior are called burn-in (or warmup). These samples don’t represent the posterior and are typically discarded.

With Metropolis-Hastings, you need to decide how many burn-in samples to discard by inspecting trace plots. NUTS handles this automatically: it uses the warmup period to adapt its step size and mass matrix, then discards those samples for you.

What NUTS gives you automatically:

  1. No tuning required: adapts step size and mass matrix during warmup.
  2. Efficient exploration: uses gradients to navigate parameter space intelligently.
  3. Automatic warmup: built-in adaptation period, no manual burn-in decisions.
  4. Better convergence: typically needs far fewer samples than Metropolis-Hastings.
TipWhich sampler should I use?
Sampler Best for Limitations
MH (Metropolis-Hastings) Simple models, discrete parameters, learning MCMC Requires manual tuning, slow for many parameters
RAM (Robust Adaptive Metropolis) Low-dimensional models where gradients aren’t available. For stochastic simulators like the particle filter we meet later in the course, gradient-free methods like this are typically what you fall back on. Still random-walk, less efficient than gradient methods
NUTS (No U-Turn Sampler) Continuous parameters with differentiable models (ODEs, most Turing models) Can’t handle discrete parameters; needs gradients

Rule of thumb: Start with NUTS(). If your model has discrete parameters or stochastic components that break automatic differentiation, fall back to MH or RAM.

Turing’s @model macro is convenient syntax, but it isn’t the only way to do MCMC in Julia. In the first session you built a log-posterior function by hand; that same function can be plugged into standalone samplers like DynamicHMC.jl via the LogDensityProblems.jl interface, giving you access to samplers beyond what Turing offers, with fewer dependencies. See Beyond Turing for a worked example.

Posterior analysis with Turing.jl

Now we can analyse the posterior distribution:

# Parameter correlations
corner_plot = marginalkde(vec(chain_nuts[:R_0]), vec(chain_nuts[:D_inf]),
                          xlabel = "R_0", ylabel = "D_inf")
┌ Warning: Keyword argument `orientation` is deprecated.
│ Please use `permute` instead.
└ @ Plots ~/.julia/packages/Plots/h49MV/src/args.jl:1558
  • marginalkde(x, y) - draws the joint density of two parameters with each marginal along its own axis, which is what a “corner plot” shows for a pair. It is how you see a correlation between parameters rather than just their separate spreads.
# Credible intervals
quantile(vec(chain_nuts[:R_0]), [0.025, 0.25, 0.5, 0.75, 0.975])
5-element Vector{Float64}:
 2.375173584839305
 2.4215239753211426
 2.444451056850713
 2.468595296358694
 2.5218742482973604

These summaries give us point estimates and uncertainty, but should we trust these results? The next session (MCMC diagnostics) teaches you how to check whether your chains converged and whether your model makes sense.

Learning points

TipLearning points
  • MCMC is a universal sampler: It works for any distribution where you can evaluate density
  • Metropolis-Hastings is simple and robust but requires tuning of proposal distributions
  • Adaptive MCMC (RAM) learns proposal covariance automatically - essential for correlated parameters
  • Turing.jl automates the mechanics of MCMC and other inference methods
  • NUTS is more efficient when gradients are available (e.g., deterministic models)
  • The Bayesian workflow: Define model → Choose sampler → Analyse samples

References

Next session: MCMC diagnostics, how to tell if your chains have converged and your results are trustworthy.

Going further (optional): MCMC is exact but slow. Variational inference is an advanced session that replaces the sampling with an optimisation. It needs only this session as background, so you can take it at any point.