Particle MCMC and model comparison

Estimated time: 75 minutes

Introduction

In the previous session, we learned how particle filters estimate likelihoods for stochastic models by simulating many trajectories and weighting them by how well they match observations. We also calibrated the number of particles needed for reliable inference.

Now we’ll put this into practice. This session shows how to combine particle filters with Metropolis-Hastings to perform full Bayesian inference on stochastic models — the Particle Marginal Metropolis-Hastings (PMMH) algorithm.

Fitting the stochastic model also puts us in a position to compare models for the first time. Two comparisons run through the second half of the session: the stochastic fit against the deterministic one, and SEITL against SEIT4L, which differ in how they distribute the duration of temporary immunity.

Coming from R?

In R, particle MCMC is available through pomp (the pmcmc() function) and the mcstate package (part of the odin/dust/mcstate ecosystem). The concepts are identical; only the implementation differs.

Objectives

In this session you will:

Implement PMMH using our calibrated particle filter from the previous session
Fit the stochastic SEIT4L model to the Tristan da Cunha outbreak
Compare stochastic vs deterministic model fits
Understand why PMMH requires long chains and adaptive proposals
Compare two competing models with DIC and with posterior predictive checks

Reminder: Stochastic variability

PMMH uses particle filters internally, introducing two sources of randomness: the MCMC proposals and the particle filter likelihood estimates. Your chains will look different from others’ runs. Focus on whether the chains have converged (diagnostics), not exact values.

Computational requirements

PMMH combines particle filters with MCMC, making it computationally demanding. For this session:

Memory: ~4 GB RAM recommended
Runtime: Each MCMC iteration takes 0.1-0.5 seconds (256 particles)
Full inference: 10,000 iterations takes ~30-60 minutes; we use pre-computed chains for the main analysis
Short test runs: the demonstration below runs 500 warmup iterations and 500 kept ones, so budget roughly twice what 500 iterations alone would cost

The pre-computed chains (loaded from data/) were generated as four chains of 162,500 iterations each, run side by side over several hours.

Setup

Source file

The source file of this session is located at sessions/pmcmc.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 FlexiChains: Parameter ## public, but not exported, so name it
using AdvancedMH ## for custom MH proposals
using CSV ## for CSV file reading
using StatsBase ## for statistical functions
using JLD2 ## for loading pre-computed chains
using SSMProblems ## for state-space model interface
using GeneralisedFilters ## for particle filtering
using ForwardDiff ## for stripping Dual numbers in non-differentiable likelihoods
using DrWatson

Load the course package

This loads the MFIIDD course package, which provides the helper functions used throughout — including the same particle filter implementation from the previous session.

using MFIIDD

SEIT4L particle filter implementation

The transitions (src/seitl_jumps.jl). Each one is a reaction with a rate, and the number of temporary immunity stages is read from the state, so the same description serves SEITL and SEIT4L:

function seitl_stoichiometry(k::Integer)
    S, E, I, L, Inc = 1, 2, 3, k + 4, k + 5
    T(j) = 3 + j  ## T_1 ... T_k sit between I and L

    reactants = [[S => 1, I => 1], [E => 1], [I => 1]]
    net = [[S => -1, E => 1], [E => -1, I => 1, Inc => 1], [I => -1, T(1) => 1]]

    for j in 1:(k - 1)  ## progress through the stages of temporary immunity
        push!(reactants, [T(j) => 1])
        push!(net, [T(j) => -1, T(j + 1) => 1])
    end

    push!(reactants, [T(k) => 1])  ## immunity wanes
    push!(net, [T(k) => -1, S => 1])

    push!(reactants, [T(k) => 1])  ## immunity becomes long term
    push!(net, [T(k) => -1, L => 1])

    return reactants, net
end

function seitl_rate_constants(θ::Dict, N::Real, k::Integer)
    β = θ[:R_0] / θ[:D_inf]
    ϵ = 1.0 / θ[:D_lat]
    ν = 1.0 / θ[:D_inf]
    τ = k / θ[:D_imm]
    α = θ[:α]

    rates = Vector{Float64}(undef, k + 4)
    rates[1], rates[2], rates[3] = β / N, ϵ, ν
    for j in 1:(k - 1)
        rates[3 + j] = τ
    end
    rates[k + 3] = (1 - α) * τ
    rates[k + 4] = α * τ
    return rates
end

SSMProblems interface (src/seitl_ssm_interface.jl):

using Random
using Distributions
using SSMProblems
using JumpProcesses: JumpProblem

"""
SEITL latent dynamics, for any number of temporary immunity stages.

State vector: `[S, E, I, T_1 ... T_k, L, daily_inc]`.

The last element carries daily incidence for the observation process. SEITL is
the `k = 1` case and SEIT4L the `k = 4` case; the number of stages comes from
the length of the state the filter is given, so one type serves both and
nothing here needs to be told which model it is running.
"""
struct SEITLDynamics{P <: JumpProblem} <: SSMProblems.LatentDynamics
    problem::P
end

"""
    SEITLDynamics(θ, init_state; rng = Random.default_rng())

Dynamics for the model whose compartments are `init_state`.

The jump problem is built once, here. Everything it needs is fixed for the whole
run: `θ`, the number of stages, and the population size, which every transition
conserves. Building it on every step doubled the cost of advancing a particle by
a day.

That makes the dynamics the owner of its random stream, so `rng` is taken here
and each `simulate` call uses it. Seed once before building the model and the
propagation is reproducible.
"""
SEITLDynamics(
    θ::Dict{Symbol, Float64},
    init_state;
    rng::AbstractRNG = Random.default_rng(),
) = SEITLDynamics(seitl_jump_problem(θ, init_state; rng = rng))

function SSMProblems.simulate(
    rng::AbstractRNG,
    dyn::SEITLDynamics,
    step::Integer,
    prev_state;
    kwargs...,
)
    ## Drop the incidence slot, leaving the compartments the simulator works on
    state = collect(prev_state[1:(end - 1)])
    return vcat(state, seitl_jump_step!(dyn.problem, state))  ## re-append incidence
end

"""
Poisson observation process.

Observes daily incidence, the last element of the state, with reporting rate ρ.
"""
struct PoissonObservation <: SSMProblems.ObservationProcess
    ρ::Float64
end

function SSMProblems.distribution(obs::PoissonObservation, step::Integer, state; kwargs...)
    daily_inc = state[end]
    return Poisson(max(obs.ρ * daily_inc, 1e-10))
end

"""
Initial state distribution, which is deterministic here: the island as the ship
lands.

Five compartments `[S, E, I, T, L]` give SEITL, eight
`[S, E, I, T1, T2, T3, T4, L]` give SEIT4L, and any other number of temporary
immunity stages works the same way.
"""
struct SEITLInitial <: SSMProblems.StatePrior
    init_state::Vector{Float64}
end

function SSMProblems.simulate(rng::AbstractRNG, prior::SEITLInitial; kwargs...)
    return vcat(prior.init_state, 0.0)  ## append 0 for initial daily incidence
end

Bootstrap filter runner (src/seitl_runner.jl):

using Base.Threads: nthreads
using Random: default_rng
using GeneralisedFilters: GeneralisedFilters, BF, DenseAncestorCallback, get_ancestry
using ForwardDiff: value
using SSMProblems: StateSpaceModel

"""
    seitl_ssm(θ, init_state)

Assemble the state-space model the two runners below both need.

`init_state` carries the number of temporary immunity stages, so this serves
SEITL (`[S, E, I, T, L]`) and SEIT4L (`[S, E, I, T1, T2, T3, T4, L]`) alike.

`init_state` is coerced to `Float64` so callers may pass an integer vector, and
ForwardDiff duals are stripped from `θ` because the bootstrap filter is not
differentiable and the component types declare concrete `Float64` fields.
"""
function seitl_ssm(θ, init_state; n_blocks::Integer = 1)
    init_f64 = collect(Float64.(init_state))
    length(init_f64) >= 5 || throw(
        ArgumentError(
            "init_state must be [S, E, I, T_1 ... T_k, L] with at least one " *
            "temporary immunity stage; got $(length(init_f64)) compartments",
        ),
    )
    θ_f64 = Dict{Symbol, Float64}(k => value(v) for (k, v) in θ)
    ## One set of dynamics per block, because a JumpProblem cannot be driven by
    ## two tasks at once. One block is the serial filter and costs nothing extra.
    dynamics = if n_blocks <= 1
        SEITLDynamics(θ_f64, init_f64)
    else
        BlockedDynamics([SEITLDynamics(θ_f64, init_f64) for _ in 1:n_blocks])
    end
    return StateSpaceModel(SEITLInitial(init_f64), dynamics, PoissonObservation(θ_f64[:ρ]))
end

"""
    run_particle_filter(θ, obs, n_particles; init_state)

Run the bootstrap particle filter and return the estimated log-likelihood.

# Arguments
- `θ`: parameter dictionary (`:R_0`, `:D_lat`, `:D_inf`, `:α`, `:D_imm`, `:ρ`)
- `obs`: vector of observed daily incidence
- `n_particles`: number of particles
- `init_state`: initial compartments; its length selects SEITL or SEIT4L
- `threaded`: propagate the particles in parallel
- `n_blocks`: how many parallel blocks to split them into

Threading pays from a few hundred particles upwards, and only once the jump
aggregation is aliased; see `seitl_jump_step!`. On six threads at 256 particles
it runs at about 1.6 times the serial speed for SEIT4L and 1.3 for SEITL, which has
fewer compartments and so less work to spread.
"""
function run_particle_filter(
    θ,
    obs,
    n_particles;
    init_state = [279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0],
    threaded::Bool = false,
    n_blocks::Integer = nthreads(),
)
    blocks = threaded ? n_blocks : 1
    model = seitl_ssm(θ, init_state; n_blocks = blocks)
    algorithm = threaded ? ThreadedBF(n_particles; n_blocks = blocks) : BF(n_particles)
    _, log_lik = GeneralisedFilters.filter(default_rng(), model, algorithm, obs)
    return log_lik
end

"""
    filtered_incidence(θ, obs, n_particles; init_state)

Run the bootstrap filter and return one draw from the smoothing distribution of
daily incidence, that is a latent path conditioned on `obs`.

`DenseAncestorCallback` records the particles and their ancestor indices at
every step, so a particle drawn from the final weights can be traced back to
give its whole path. The last element of the state is the daily incidence, so
reading that element off the path gives the trajectory directly.

# Returns
- `Vector{Float64}` of length `length(obs)`, the filtered daily incidence
"""
function filtered_incidence(
    θ,
    obs,
    n_particles;
    init_state = [279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0],
)
    model = seitl_ssm(θ, init_state)
    callback = DenseAncestorCallback(nothing)
    final, _ =
        GeneralisedFilters.filter(default_rng(), model, BF(n_particles), obs; callback)

    ## draw one particle in proportion to its final weight, then follow its
    ## ancestry back to the start
    log_w = getfield.(final.particles, :log_w)
    w = exp.(log_w .- maximum(log_w))
    u = rand() * sum(w)
    idx = findfirst(>=(u), cumsum(w))
    path = get_ancestry(callback.container, idx)

    return [path[t][end] for t in 1:length(obs)]
end

Initialisation

We set a random seed for reproducibility.

Random.seed!(1234)

TaskLocalRNG()

Recap: Data and particle filter

We load the Tristan da Cunha data and set up initial conditions:

# Load Tristan da Cunha data
flu_tdc = CSV.read(datadir("flu_tdc_1971.csv"), DataFrame)

# Initial state for SEIT4L
init_state = Dict(
    :S => 279.0,
    :E => 0.0,
    :I => 2.0,
    :T1 => 3.0,
    :T2 => 0.0,
    :T3 => 0.0,
    :T4 => 0.0,
    :L => 0.0
)

Dict{Symbol, Float64} with 8 entries:
  :I  => 2.0
  :T2 => 0.0
  :L  => 0.0
  :S  => 279.0
  :T4 => 0.0
  :T1 => 3.0
  :E  => 0.0
  :T3 => 0.0

In the Particle filters session, we:

Wrote the SEITL transitions once, for any number of temporary immunity stages (seitl_stoichiometry, seitl_rate_constants)
Defined the SSMProblems.jl interface (SEITLDynamics, PoissonObservation, SEITLInitial)
Calibrated the particle count against the SD of the log-likelihood estimate, which put the minimum acceptable count at 64

All these components were loaded via using MFIIDD in the setup above (defined in src/seitl_jumps.jl, src/seitl_ssm_interface.jl, and src/seitl_runner.jl). Now we’ll use them for full Bayesian inference.

In your own work, you should run prior and posterior predictive checks when fitting models. We skip them for now to focus on the PMMH algorithm.

Particle MCMC (PMMH) with Turing.jl

Now that we have a calibrated particle filter, we can use it for inference. We’ll implement Particle Marginal Metropolis-Hastings (PMMH) - one of the two main particle MCMC algorithms from Andrieu et al. (2010).

PMMH uses the particle filter to estimate the likelihood, then does standard Metropolis-Hastings on the parameters. In Turing.jl, we implement this by:

Defining priors as usual with ~
Building a state-space model and running GeneralisedFilters.filter() to get the log-likelihood
Adding it to the model’s log-probability with @addlogprob!
Sampling with MH() with custom random-walk proposals

"""
PMMH model for SEIT4L inference.
Uses GeneralisedFilters for likelihood estimation.
"""
@model function pmmh_seit4l(obs, n_particles)
    # Weakly informative priors
    R_0 ~ truncated(Normal(3.0, 2.0), lower=1.0)
    D_lat ~ truncated(Normal(2.0, 1.0), lower=0.5)
    D_inf ~ truncated(Normal(3.0, 2.0), lower=0.5)
    α ~ Beta(2, 2)
    D_imm ~ truncated(Normal(15.0, 10.0), lower=1.0)
    ρ ~ Beta(2, 2)

    # Create parameter dictionary; extract underlying values so ForwardDiff
    # Duals (used only by Turing's gradient-based init probe) don't flow into
    # the particle filter, which is non-differentiable anyway.
    θ = Dict(:R_0 => ForwardDiff.value(R_0),
             :D_lat => ForwardDiff.value(D_lat),
             :D_inf => ForwardDiff.value(D_inf),
             :α => ForwardDiff.value(α),
             :D_imm => ForwardDiff.value(D_imm),
             :ρ => ForwardDiff.value(ρ))

    # Build state-space model
    # The SSMProblems interface expects a plain vector, not a Dict
    # Order: [S, E, I, T1, T2, T3, T4, L]
    init_state = [279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0]
    ssm = StateSpaceModel(SEITLInitial(init_state),
                          SEITLDynamics(θ, init_state),
                          PoissonObservation(θ[:ρ]))

    # Run particle filter to estimate log-likelihood
    _, log_lik = GeneralisedFilters.filter(Random.default_rng(), ssm, BF(n_particles), obs)

    # Add log-likelihood to model (priors are already included via ~)
    Turing.@addlogprob! log_lik
end

Main.Notebook.pmmh_seit4l

Code explainer

StateSpaceModel(...) - assembles the state-space model from initial distribution, dynamics, and observation components defined in src/seitl_ssm_interface.jl.
GeneralisedFilters.filter(...) - runs the bootstrap particle filter (BF) and returns the estimated log-likelihood. Random.default_rng() provides the global random number generator — this is Julia’s way of passing randomness explicitly rather than relying on a hidden global state.
@addlogprob! - manually adds the particle filter’s log-likelihood estimate to the model. The priors specified with ~ still contribute automatically. This is the recommended pattern for integrating external likelihood estimators with Turing.
Why 256 particles? PMMH acceptance depends on the variance of the log-likelihood estimate. Using more particles reduces variance but increases computation. The calibration in the previous session put the minimum at 64, measured at a single parameter set; how noisy the estimate is depends on where in the posterior you are, so it is safer to go somewhat higher, and worth checking the count again once you have a posterior to check it against.

Now let’s run PMMH. Unlike standard MCMC where NUTS auto-tunes proposals, PMMH requires Metropolis-Hastings because the particle filter likelihood is not differentiable.

Let’s start with a short chain to see the challenges:

# Define model with 256 particles
model_pmmh = pmmh_seit4l(flu_tdc.obs, 256)

# Run a short chain with RAM
# externalsampler() wraps an AdvancedMH.jl sampler for use with Turing models
# num_warmup=500 gives it 500 discarded iterations to adapt in
# check_model=false is needed because our model uses @addlogprob! instead of standard ~ syntax
chain_short = sample(
    model_pmmh,
    externalsampler(AdvancedMH.RobustAdaptiveMetropolis()),
    500;
    num_warmup=500,
    check_model=false,
    progress=false
)

╭─FlexiChain (500 iterations, 1 chain) ────────────────────────────────────────╮
│ ↓ iter  = 501:1000                                                           │
│ → chain = 1:1                                                                │
│                                                                              │
│ Parameters (6) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_lat, D_inf, α, D_imm, ρ                                     │
│                                                                              │
│ Extras (6)                                                                   │
│  Float64  logα, η, logprior, loglikelihood, logjoint                         │
│  Bool     accepted                                                           │
╰──────────────────────────────────────────────────────────────────────────────╯

# Plot the short chain
plot(chain_short)

Notice the poor mixing

Look carefully at the trace plots above. You should see:

Flat stretches (plateaus): The chain gets stuck at the same value for many iterations
Slow wandering: Instead of rapidly exploring, the chain drifts gradually

With only 500 samples, the chain hasn’t had time to properly explore the posterior.

Two details of the figure are worth reading. The iteration axis starts at 500 because the 500 warmup iterations are discarded, so what you are looking at is the chain after RAM has had a go at adapting. And it is still stuck: adaptation needs far more than 500 iterations to find a proposal that suits this posterior, so those flat stretches are not simply a badly chosen step size that a little tuning would cure.

Exercise: Run your own short PMMH chain

Try running a short PMMH chain yourself:

Run sample(model_pmmh, externalsampler(AdvancedMH.RobustAdaptiveMetropolis()), 200; num_warmup=200, check_model=false) and look at the trace plots
What acceptance rate does the sampler achieve? (Look at the proportion of samples that change between iterations)
Try reducing the particle count to 64: pmmh_seit4l(flu_tdc.obs, 64). How does this affect mixing? Is it faster but noisier?
Compare to the pre-computed chain loaded below — how many more iterations were needed for convergence?

What to look for

Expect a low acceptance rate and a trace plot with long flat stretches — the chain sitting at one parameter value for many iterations before moving. That is the signature of a noisy likelihood: once the sampler lands on a point whose likelihood was over-estimated by chance, every honest proposal has to beat an inflated number, so it stays put.

Dropping to 64 particles makes each iteration roughly twice as fast and the mixing distinctly worse, because halving the particles increases the variance of the estimate. Whether that is a good trade depends on which side of the optimum you started, which is the calibration argument from the previous session in practical form.

Against the pre-computed chain, 200 iterations is nowhere near enough. The saved chain runs orders of magnitude longer, and comparing your trace with it makes the point better than any convergence statistic: a short PMMH chain does not look obviously broken, it just has not gone anywhere yet.

Why is PMMH mixing difficult?

The particle filter gives a noisy estimate of the likelihood. Even at the same parameters, different runs give different values. This noise means proposals that would be accepted with an exact likelihood may be rejected (and vice versa), slowing exploration.

Robust Adaptive Metropolis (RAM) helps by learning the posterior covariance structure during its warmup iterations (Vihola, 2012), and then sampling with the proposal it has learned.

Pre-computed chains

Running PMMH with RAM for good mixing takes several hours. We load pre-computed chains here: four chains run side by side, each with 50,000 warmup iterations for RAM to adapt in, then 112,500 kept and thinned by 50, which is 9,000 draws in total. Four chains rather than one long chain because they run in the time of one, and because R-hat compares chains that started independently. See scripts/generate_pmcmc_seit4l.jl if you want to run it yourself.

# Load the pre-computed SEIT4L chains with good mixing
chain_df = CSV.read(datadir("pmcmc_seit4l_chain.csv"), DataFrame)
chain_params = [:R_0, :D_lat, :D_inf, :α, :D_imm, :ρ]

# FlexiChains builds a chain from an iterations × chains × parameters array, keeping
# the parameters in the order given. The file holds the four chains stacked, one
# per value of `chain`, so splitting on that column and stacking along the second
# dimension is what lets R-hat compare separate chains
draws = stack([Matrix(g[:, chain_params]) for g in groupby(chain_df, :chain)]; dims = 2)
chain_pmcmc = SymChain(draws, Tuple(Parameter.(chain_params)))

╭─FlexiChain (2250 iterations, 4 chains) ──────────────────────────────────────╮
│ ↓ iter  = 1:2250                                                             │
│ → chain = 1:4                                                                │
│                                                                              │
│ Parameters (6) ── Symbol                                                     │
│  Float64  R_0, D_lat, D_inf, α, D_imm, ρ                                     │
│                                                                              │
│ Extras (0)                                                                   │
│  (none)                                                                      │
╰──────────────────────────────────────────────────────────────────────────────╯

# Display chain summary
summarystats(chain_pmcmc)

╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (6) ── Symbol                                                     │
│  Float64  R_0, D_lat, D_inf, α, D_imm, ρ                                     │
│                                                                              │
│ Extras (0)                                                                   │
│  (none)                                                                      │
│                                                                              │
│ Summary                                                                      │
│   param     mean     std    mcse   ess_bulk   ess_tail    rhat      q5  …    │
│     R_0   6.3176  1.1620  0.0155  5555.5178  7353.6998  1.0008  4.5310  …    │
│   D_lat   1.3365  0.4167  0.0091  1794.9129  1124.0896  0.9999  0.6565  …    │
│   D_inf   2.1151  0.8189  0.0165  2186.9533  1809.2534  1.0002  0.8153  …    │
│       α   0.4703  0.0553  0.0006  8121.8009  7609.8112  0.9999  0.3756  …    │
│   D_imm  11.7624  1.8503  0.0207  7965.2780  7498.6460  1.0001  8.8557  …    │
│       ρ   0.6921  0.0481  0.0005  8399.4977  8622.4198  1.0004  0.6138  …    │
╰──────────────────────────────────────────────────────────────────────────────╯

# Plot trace plots - much better mixing!
plot(chain_pmcmc)

Compare these traces to the short chain above. The chains cross the posterior repeatedly instead of sticking for hundreds of iterations, and the summary above puts numbers on it. Read R-hat first: it compares the variance within each of the four chains against the variance between them, so a value at 1.00 says four chains that started in different places agree about where the posterior is. The effective sample size then says how many independent draws the 9,000 kept ones are worth, which for a correlated chain is always fewer: here the worst parameter is $D_\mathrm{lat}$ at about 1,800, and $\rho$ , which the data pin down hardest, is worth over 8,000.

The saved chain also records whether each draw was accepted, which is the diagnostic that says whether RAM’s adaptation worked:

# `accepted` is kept alongside the parameters, so we can check the rate RAM reached
println("acceptance rate: ", round(100 * mean(chain_df.accepted), digits = 1), "%")

acceptance rate: 22.9%

That is the target RAM aims for, and it is the payoff from the warmup iterations described below.

The case for robust adaptive Metropolis

With standard MH and independent proposals, mixing is poor when parameters are correlated. The sampler proposes moves that go “against” the correlation structure and gets rejected.

RAM solves this by learning the posterior covariance (Vihola, 2012). During warmup it:

Adapts the proposal covariance to match the posterior shape
Scales proposals towards ~23% acceptance (optimal for random-walk MH)

The adaptation happens in the warmup iterations only, with a step size that decays as $n^{-0.6}$ ; at the end of warmup the proposal is fixed and the chain that follows is an ordinary Metropolis chain. This is similar to how NUTS adapts its mass matrix during warmup, but works for non-differentiable likelihoods like particle filters.

Warmup is the whole adaptation budget

Because RAM only adapts during warmup, num_warmup is not a formality: with no warmup the proposal stays at the identity matrix it started from and acceptance sits near 0.5% instead of the 23% RAM aims for, whatever the chain length. The saved chains use 50,000 warmup iterations.

Comparing with deterministic model fits

We quickly fit the deterministic SEIT4L model for comparison:

"""
SEIT4L deterministic model.
"""
function seit4l_ode!(du, u, p, t)
    R_0, D_lat, D_inf, α, D_imm = p
    β = R_0 / D_inf
    ϵ = 1.0 / D_lat
    ν = 1.0 / D_inf
    τ = 4.0 / D_imm

    S, E, I, T1, T2, T3, T4, L, Inc = u
    N = S + E + I + T1 + T2 + T3 + T4 + L

    du[1] = -β * S * I / N + (1 - α) * τ * T4
    du[2] = β * S * I / N - ϵ * E
    du[3] = ϵ * E - ν * I
    du[4] = ν * I - τ * T1
    du[5] = τ * T1 - τ * T2
    du[6] = τ * T2 - τ * T3
    du[7] = τ * T3 - τ * T4
    du[8] = α * τ * T4
    du[9] = ϵ * E
end

seit4l_prob = ODEProblem(seit4l_ode!, zeros(9), (0.0, 1.0), zeros(5))

@model function seit4l_deterministic_model(init_state, times, n_obs)
    # Weakly informative priors
    R_0 ~ truncated(Normal(3.0, 2.0), lower=1.0)
    D_lat ~ truncated(Normal(2.0, 1.0), lower=0.5)
    D_inf ~ truncated(Normal(3.0, 2.0), lower=0.5)
    α ~ Beta(2, 2)
    D_imm ~ truncated(Normal(15.0, 10.0), lower=1.0)
    ρ ~ Beta(2, 2)

    # Solve ODE; promote u0 to the parameter element type so ForwardDiff Duals
    # propagate through the solution
    params = [R_0, D_lat, D_inf, α, D_imm]
    T = eltype(params)
    u0 = T[init_state[:S], init_state[:E], init_state[:I],
           init_state[:T1], init_state[:T2], init_state[:T3], init_state[:T4],
           init_state[:L], 0.0]

    prob = remake(
        seit4l_prob, u0 = u0, p = params, tspan = (times[1], times[end])
    )
    sol = solve(prob, Tsit5(), saveat=times)

    # Reject parameter draws that drive the solver to fail early; the gradient
    # probe in NUTS init occasionally lands on combinations where the ODE only
    # saves the initial timepoint.
    if !SciMLBase.successful_retcode(sol) || length(sol.t) != length(times)
        Turing.@addlogprob! -Inf
        return
    end

    # Extract daily incidence. There is one fewer difference than there are time
    # points, and `inc_daily[i]` is the incidence over the day ending at
    # `times[i + 1]`, which is the day the i-th observation counts. Prepending a
    # zero instead would match every observation to the previous day.
    inc_cumulative = sol[9, :]
    inc_daily = diff(inc_cumulative)

    # Likelihood: Poisson observations
    lambdas = [max(ρ * inc_daily[i], 1e-10) for i in 1:n_obs]
    obs ~ arraydist(Poisson.(lambdas))
end

seit4l_deterministic_model (generic function with 2 methods)

Why | conditioning here but @addlogprob! above?

The deterministic model uses obs ~ arraydist(Poisson.(lambdas)) — the likelihood is expressed as a ~ statement, so Turing can condition on the data directly with |.

The stochastic pmmh_seit4l model uses an external likelihood estimator (the particle filter via GeneralisedFilters.filter), which returns a log-likelihood that we add manually with @addlogprob!. In that model, there is no obs ~ statement in the model body — obs enters as a function argument passed to the filter — so | has nothing to condition on.

In short: use | when the likelihood is a ~ statement; use @addlogprob! when the likelihood comes from an external computation.

# Fit deterministic model
times_obs = 0.0:1.0:59.0
model_det = seit4l_deterministic_model(init_state, times_obs, length(flu_tdc.obs)) | (; obs = flu_tdc.obs)
chain_det = sample(model_det, NUTS(), 2000, progress=false)

┌ Info: Found initial step size
└   ϵ = 0.2
┌ Warning: Verbosity toggle: dt_epsilon 
│  Initial timestep too small (near machine epsilon), using default: dt = 1.0e-6
└ @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/luSGN/src/initdt.jl:222
┌ Warning: Verbosity toggle: max_iters 
│  Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/iMnPY/src/check_error.jl:20

╭─FlexiChain (2000 iterations, 1 chain) ───────────────────────────────────────╮
│ ↓ iter  = 1001:3000                                                          │
│ → chain = 1:1                                                                │
│                                                                              │
│ Parameters (6) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_lat, D_inf, α, D_imm, ρ                                     │
│                                                                              │
│ Extras (14)                                                                  │
│  Int64    n_steps, tree_depth                                                │
│  Bool     is_accept, numerical_error                                         │
│  Float64  acceptance_rate, log_density, hamiltonian_energy,                  │
│           hamiltonian_energy_error, max_hamiltonian_energy_error, step_size, │
│           nom_step_size, logprior, loglikelihood, logjoint                   │
╰──────────────────────────────────────────────────────────────────────────────╯

Plotting code

## Compare posteriors, both models overlaid on one panel per parameter
params = [:R_0, :D_lat, :D_inf, :α, :D_imm, :ρ]
panels = map(enumerate(params)) do (i, par)
    ## pull the draws out as plain vectors so the series labels are ours
    pl = density(vec(chain_det[par]), title = string(par), label = "deterministic",
                 colour = :steelblue, legend = i == 1 ? :topright : false)
    density!(pl, vec(chain_pmcmc[par]), label = "stochastic", colour = :darkorange)
end

plot(panels..., layout = (2, 3), size = (900, 500))

Exercise: Interpret the stochastic vs deterministic comparison

Look at the posterior density plots above and compare the two models:

Which model gives wider posterior distributions? Why would accounting for stochastic variability change the uncertainty in parameter estimates?
Focus on the $\alpha$ parameter (probability of long-term immunity). Do the two models agree on its value?
For a population of 284 people, do you think the deterministic approximation is reasonable, or does stochasticity matter?

What to look for

The stochastic model gives wider posteriors, but by very different amounts across the six parameters. The clearest case is $D_\mathrm{imm}$ , whose posterior standard deviation roughly doubles, followed by $\alpha$ and $D_\mathrm{lat}$ , which grow by about three fifths. $R_0$ grows by about a tenth, and the reporting rate $\rho$ barely moves at all, which makes sense: it is pinned by the total number of cases either way.

The direction is the same everywhere, and it should be. The deterministic model attributes every discrepancy between prediction and data to observation noise, so the parameters have to explain everything and get pinned down tightly. The stochastic model has a second source of variation available — the epidemic itself could have gone differently — so a wider range of parameters remains consistent with what was seen. A narrower posterior is not a better one when the narrowness comes from a mechanism you left out.

On $\alpha$ , look at whether the intervals overlap rather than whether the means match. Disagreement here matters more than for the other parameters, because $\alpha$ controls how much of the population returns to susceptible and therefore whether a second wave happens at all.

On 284 people: stochasticity matters. Demographic noise scales roughly as $1/\sqrt{N}$ , so a few hundred people is squarely in the regime where chance events change the trajectory — and Tristan da Cunha is the standard example precisely because the population is small enough for it to show. The deterministic approximation earns its keep in populations of hundreds of thousands, not hundreds.

For the Tristan da Cunha outbreak with only 284 individuals, demographic stochasticity plays a significant role. The stochastic model typically shows wider posteriors because it accounts for an additional source of variability — the randomness of individual infection events — that the deterministic model ignores. The deterministic model is much faster to fit and useful for initial exploration, but for small populations the stochastic model gives more honest uncertainty estimates.

Model comparison: SEITL vs SEIT4L

In the Structure and stochasticity session, we introduced two model variants that both explain the two-wave pattern through temporary immunity:

SEITL: Single T compartment with exponential duration (memoryless)
SEIT4L: Four T sub-compartments giving Erlang-4 duration (more realistic timing)

Which model better explains the data? We fit both and compare.

Fitting SEITL for comparison

The SEITL model uses the same structure as SEIT4L, but with a single T compartment instead of four. The implementation follows the same pattern we used for SEIT4L: a jump simulator for stochastic dynamics, a bootstrap particle filter to estimate the likelihood, and a Turing @model that adds the log-likelihood via @addlogprob!. Since the code structure is identical, we load pre-computed chains rather than repeating the implementation. See scripts/generate_pmcmc_seitl.jl if you want to examine the full code.

Both models were fitted using RAM as four chains each, with 50,000 warmup iterations per chain, then 112,500 kept per chain and thinned by 50:

# Load pre-computed SEITL chain
chain_seitl_df = CSV.read(datadir("pmcmc_seitl_chain.csv"), DataFrame)
draws_seitl =
    stack([Matrix(g[:, chain_params]) for g in groupby(chain_seitl_df, :chain)]; dims = 2)
chain_seitl = SymChain(draws_seitl, Tuple(Parameter.(chain_params)))

╭─FlexiChain (2250 iterations, 4 chains) ──────────────────────────────────────╮
│ ↓ iter  = 1:2250                                                             │
│ → chain = 1:4                                                                │
│                                                                              │
│ Parameters (6) ── Symbol                                                     │
│  Float64  R_0, D_lat, D_inf, α, D_imm, ρ                                     │
│                                                                              │
│ Extras (0)                                                                   │
│  (none)                                                                      │
╰──────────────────────────────────────────────────────────────────────────────╯

# Compare the fits visually
p1 = density(chain_seitl, :R_0, title="R_0", label="SEITL")
density!(p1, chain_pmcmc, :R_0, label="SEIT4L")

p2 = density(chain_seitl, :D_imm, title="D_imm", label="SEITL")
density!(p2, chain_pmcmc, :D_imm, label="SEIT4L")

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

Approaches to model comparison

Looking at the posterior densities above, we can see differences between the two models. But which model better explains the data? This is the question of model comparison or model selection.

There are three families of answer, and which of them applies is settled by how the model is built.

Information criteria estimate how well a model would predict new data, with a penalty for complexity. The ones usually recommended, the widely applicable information criterion (WAIC) and leave-one-out cross-validation (LOO), need a likelihood term for each observation on its own, and a stochastic model does not give you those: every day’s count depends on the same random trajectory, so a day can only be scored given the days before it. The Deviance Information Criterion (DIC) needs only the likelihood of the whole series, which is what the particle filter estimates. However, it is a rough approximation to predictive accuracy, so do interpret its results with caution.

Bayes factors compare marginal likelihoods directly:

$\text{BF}_{12} = \frac{P(\text{data}|M_1)}{P(\text{data}|M_2)}$

Values > 10 indicate strong evidence for $M_1$ ; values < 0.1 indicate strong evidence for $M_2$ . The marginal likelihood can be estimated from particle filter output, for instance by averaging the filter’s likelihood estimates over draws from the prior, but this is computationally intensive.

Posterior predictive checks ask which model generates data that looks like the observations. They give no single number, and what they offer instead is a view of how the models differ. They also apply whatever the model’s structure, which is the reason this session leans on them.

Comparing parameter estimates

We can compare how the two models estimate the same parameters:

# Posterior means from each model
println("Posterior mean estimates:")
println("SEITL:")
println("  R_0 = $(round(mean(chain_seitl[:R_0]), digits=2))")
println("  D_imm = $(round(mean(chain_seitl[:D_imm]), digits=2))")

println("SEIT4L:")
println("  R_0 = $(round(mean(chain_pmcmc[:R_0]), digits=2))")
println("  D_imm = $(round(mean(chain_pmcmc[:D_imm]), digits=2))")

Posterior mean estimates:
SEITL:
  R_0 = 5.72
  D_imm = 15.29
SEIT4L:
  R_0 = 6.32
  D_imm = 11.76

The key difference is in the D_imm (duration of temporary immunity) estimates. SEIT4L, with its Erlang-distributed immunity duration, typically estimates a shorter mean duration but with less variance in timing.

Comparing the two models

The criteria above are described but not computed, so nothing so far has actually compared the two models. We can do it directly with the chains we already have.

DIC needs the deviance $D(\theta) = -2 \log p(y \mid \theta)$ , averaged over the posterior and at the posterior mean, and both come from the likelihood of the whole series. That is what a particle filter gives you, and it is why Camacho et al. (2011) used DIC on these models.

First, the log-likelihood of every draw is already saved:

seitl_ll  = CSV.read(datadir("pmcmc_seitl_chain.csv"), DataFrame).loglikelihood
seit4l_ll = CSV.read(datadir("pmcmc_seit4l_chain.csv"), DataFrame).loglikelihood
println("mean log-likelihood, SEITL:  ", round(mean(seitl_ll), digits = 2))
println("mean log-likelihood, SEIT4L: ", round(mean(seit4l_ll), digits = 2))

mean log-likelihood, SEITL:  -120.94
mean log-likelihood, SEIT4L: -117.77

SEIT4L is ahead by about three log-likelihood units, more than the spread of either chain, so this is not sampling noise. That agrees with Camacho et al. (2011), who preferred SEIT4L.

On its own that is not a comparison: a mean log-likelihood rewards fit and penalises nothing, so it will always favour the more flexible model. DIC adds a penalty for complexity. Writing $\overline{D}$ for the deviance averaged over the posterior and $\bar{\theta}$ for the posterior mean,

$\mathrm{DIC} = D(\bar{\theta}) + 2 p_D, \qquad p_D = \overline{D} - D(\bar{\theta}),$

which is the same as $\overline{D} + p_D$ . The three quantities are:

$D(\bar{\theta})$ , the fit at the posterior mean, a single best guess at the parameters.
$\overline{D}$ , the fit averaged over everything the posterior still considers plausible. It is minus twice the mean log-likelihood over the posterior, re-estimated below rather than read off the draws, for the reason given just before the cell.
$p_D$ , the effective number of parameters: how much worse the average fit is than the fit at the mean. Each parameter the data pin down adds about one, and a parameter the data leave to its prior adds close to nothing, so $p_D$ need not be a whole number and falls below the number of parameters you wrote down.

Both deviances need fresh runs through the particle filter of the model that produced the chain. The saved log-likelihood is not the one to average for $\overline{D}$ : it is the estimate that was current when the draw was recorded, and a draw is accepted partly because its estimate came out high, so the saved values sit above a fresh run by about half a log unit. Re-estimating over a sample of draws puts $\overline{D}$ and $D(\bar{\theta})$ on the same footing, which is what makes the filter’s own bias cancel in $p_D$ .

## DIC from a chain and the particle filter of the model it was fitted with
const DIC_DRAWS = 200  ## D̄ to about a third of a deviance unit

function dic(chain_df, run_filter)
    params = (:R_0, :D_lat, :D_inf, :α, :D_imm, :ρ)

    ## D̄: one fresh run at each of a sample of draws
    rows = rand(1:nrow(chain_df), DIC_DRAWS)
    ll_post = [
        run_filter(Dict(p => chain_df[i, p] for p in params), flu_tdc.obs, 256) for i in rows
    ]
    D_mean = -2 * mean(ll_post)

    θ_mean = Dict(p => mean(chain_df[!, p]) for p in params)
    ## the filter estimate is noisy, so average several runs at the posterior mean
    ll = [run_filter(θ_mean, flu_tdc.obs, 256) for _ in 1:20]
    D_at_mean = -2 * mean(ll)
    p_D = D_mean - D_at_mean
    return (; DIC = D_at_mean + 2 * p_D, D_mean, D_at_mean, p_D, var_ll = var(ll))
end

Random.seed!(1)
## One filter serves both. Five compartments give SEITL, eight give SEIT4L.
for (name, file, init) in (
    ("SEITL", "pmcmc_seitl_chain.csv", [279.0, 0.0, 2.0, 3.0, 0.0]),
    ("SEIT4L", "pmcmc_seit4l_chain.csv", [279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0]),
)
    run_filter = (θ, o, n) -> run_particle_filter(θ, o, n; init_state = init)
    d = dic(CSV.read(datadir(file), DataFrame), run_filter)
    println(rpad(name, 8), "D̄ = ", round(d.D_mean, digits = 1),
            "  D(θ̄) = ", round(d.D_at_mean, digits = 1),
            "  p_D = ", round(d.p_D, digits = 1),
            "  DIC = ", round(d.DIC, digits = 1),
            "  var(log p̂) at θ̄ = ", round(d.var_ll, digits = 2))
end

SEITL   D̄ = 242.8  D(θ̄) = 238.5  p_D = 4.3  DIC = 247.2  var(log p̂) at θ̄ = 0.46
SEIT4L  D̄ = 235.9  D(θ̄) = 231.3  p_D = 4.6  DIC = 240.5  var(log p̂) at θ̄ = 0.34

Lower is better, and only the difference between two models means anything. SEIT4L comes out about seven lower, and all of the difference is fit. The two models have about the same effective number of parameters, which is what you would expect: they have the same six parameters, and SEIT4L’s four immunity compartments share a single duration, so they add none. What the four compartments change is the shape of the time spent immune, and that is what the data reward.

There are two caveats that apply to using DIC here. Both deviances come from particle filter estimates, and a filter run underestimates $\log p(y \mid \theta)$ on average, so each deviance comes out a little high. That applies to $\overline{D}$ as much as to $D(\bar{\theta})$ , because the cell above re-estimates both with fresh runs of the same filter at the same particle count. The bias therefore largely cancels in $p_D$ , which is why it is worth paying for those extra runs rather than averaging the log-likelihood stored with each draw. The noise does not: a noisier estimate at the posterior mean moves DIC directly, which is why the cell averages twenty runs there, and the variance it prints shows how far a single run wanders.

The larger one is that DIC is not uniquely defined for a model with latent variables (Celeux et al. 2006). It can be built on the likelihood with the latent path integrated out, or with the path treated as a parameter, and the two give different answers. The version above is the first, which is what the particle filter estimates. The second does not suit a model like this one: it carries one latent state per observation, so its effective number of parameters grows with the data. Other ways to compare these models are listed under Going further.

No single number can say whether the winner fits the data at all, so the comparison this session leans on is the posterior predictive one. Simulate from each posterior and look:

Each replicate draws parameters from the posterior and then runs the epidemic forward from the same starting state, so it gets a fresh latent trajectory. That costs wide envelopes, since every replicate re-rolls the whole outbreak along with the parameters. The alternative is to condition the trajectory on the observations, and we compare the two below.

init_seitl = Dict(:S => 279.0, :E => 0.0, :I => 2.0, :T => 3.0, :L => 0.0)
init_4l = [279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0]
n_obs = length(flu_tdc.obs)

## one stochastic replicate per posterior draw, observed through the same
## Poisson reporting model the filter used
function predict_seitl(row)
    θ = Dict(:R_0 => row.R_0, :D_lat => row.D_lat, :D_inf => row.D_inf,
             :α => row.α, :D_imm => row.D_imm)
    traj = simulate_seitl_stochastic(θ, init_seitl, times_obs)
    [rand(Poisson(max(row.ρ * traj.Inc[i + 1], 1e-10))) for i in 1:n_obs]
end

function predict_seit4l(row)
    θ = Dict(:R_0 => row.R_0, :D_lat => row.D_lat, :D_inf => row.D_inf,
             :α => row.α, :D_imm => row.D_imm)
    state = copy(init_4l)
    problem = seitl_jump_problem(θ, state)
    [rand(Poisson(max(row.ρ * seitl_jump_step!(problem, state, 1.0), 1e-10)))
     for _ in 1:n_obs]
end

Random.seed!(1)
n_sims = 200
seitl_df = CSV.read(datadir("pmcmc_seitl_chain.csv"), DataFrame)
seit4l_df = CSV.read(datadir("pmcmc_seit4l_chain.csv"), DataFrame)
sims_seitl = [predict_seitl(seitl_df[rand(1:nrow(seitl_df)), :]) for _ in 1:n_sims]
sims_seit4l = [predict_seit4l(seit4l_df[rand(1:nrow(seit4l_df)), :]) for _ in 1:n_sims]

function envelope!(p, sims, label, colour)
    M = reduce(hcat, sims)
    lo = [quantile(M[t, :], 0.025) for t in 1:size(M, 1)]
    hi = [quantile(M[t, :], 0.975) for t in 1:size(M, 1)]
    med = [median(M[t, :]) for t in 1:size(M, 1)]
    plot!(p, flu_tdc.time, med, ribbon = (med .- lo, hi .- med),
          label = label, colour = colour, fillalpha = 0.2, linewidth = 2)
end

p = plot(xlabel = "Day", ylabel = "Reported cases",
         title = "Posterior predictive, 95% envelopes")
envelope!(p, sims_seitl, "SEITL", :steelblue)
envelope!(p, sims_seit4l, "SEIT4L", :darkorange)
scatter!(p, flu_tdc.time, flu_tdc.obs, label = "observed",
         colour = :black, markersize = 3)

Both envelopes contain almost every observation, so neither model is refuted by the data as a whole. The interesting question is the second wave, which the two models disagree about mechanically. Take its size as a test statistic and compute a Bayesian p-value, as in the model checking session:

second_wave = 25:45
observed_second = sum(flu_tdc.obs[second_wave])
for (name, sims) in (("SEITL", sims_seitl), ("SEIT4L", sims_seit4l))
    totals = [sum(v[second_wave]) for v in sims]
    println(rpad(name, 7), " second-wave cases: median ", median(totals),
            ", p = ", round(mean(totals .>= observed_second), digits = 3))
end
println("observed: ", observed_second)

SEITL   second-wave cases: median 64.0, p = 0.13
SEIT4L  second-wave cases: median 69.0, p = 0.225
observed: 88

What to look for

SEIT4L generates the larger second wave, which is the direction its mechanism predicts: Erlang-distributed immunity returns people to susceptible in a reasonably synchronised group, where SEITL’s exponential waiting time returns them steadily from the moment they recover and smears the wave out.

Notice what the p-values do not say, though. Both models undershoot the observed second wave, and neither is extreme enough to be rejected on this statistic. The better model here is the less-bad one. That is the limit of a comparison: it ranks the candidates you happened to write down, and says nothing about whether any of them is right.

Where the latent trajectory comes from

The replicates above simulate the epidemic forward, so each one carries a latent path the filter never saw the data for. A particle filter can also hand back a path that was conditioned on the data. DenseAncestorCallback records the particles and their ancestor indices at every step, so drawing one particle at the end and following its ancestry back gives a whole trajectory:

Random.seed!(2)
n_filt = 200
sims_filtered = map(1:n_filt) do _
    row = seit4l_df[rand(1:nrow(seit4l_df)), :]
    θ = Dict(:R_0 => row.R_0, :D_lat => row.D_lat, :D_inf => row.D_inf,
             :α => row.α, :D_imm => row.D_imm, :ρ => row.ρ)
    inc = filtered_incidence(θ, flu_tdc.obs, 256)
    [rand(Poisson(max(row.ρ * x, 1e-10))) for x in inc]
end

p = plot(xlabel = "Day", ylabel = "Reported cases",
         title = "SEIT4L: re-simulated and filtered, 95% envelopes")
envelope!(p, sims_seit4l, "re-simulated", :darkorange)
envelope!(p, sims_filtered, "filtered", :seagreen)
scatter!(p, flu_tdc.time, flu_tdc.obs, label = "observed",
         colour = :black, markersize = 3)

Both sets are draws of replicate data from the same fitted model, and they differ only in where the latent path comes from. The re-simulated replicates draw $\theta \sim p(\theta \mid y)$ and then a fresh path $X \sim p(X \mid \theta)$ , with the latent path marginalised out of what the data informed. The filtered replicates draw the pair from the joint $p(\theta, X \mid y)$ , so the path itself has been told the answer.

Warning

The filtered envelope is narrower on 51 of the 59 days, and less than half the width of the re-simulated one across the first wave. That tightness is not evidence of a better model. It tracks the observations because it was conditioned on them, and it would do so for a model that generates nothing like this outbreak when left to run. Judge fit with the re-simulated envelope. Read the filtered one as a picture of what the epidemic plausibly did, which is the question it answers well.

Findings of the original study

Camacho et al. (2011) compared SEITL and SEIT4L using DIC and found that SEIT4L provided a better fit (lower DIC). This suggests that the timing of immunity loss matters - it’s not memoryless, as the exponential distribution in SEITL would imply.

The Erlang-4 distribution in SEIT4L has a coefficient of variation of 0.5 (compared to 1.0 for exponential), meaning individuals lose immunity at more consistent times. For the Tristan da Cunha data, this more realistic timing distribution better captures the second wave dynamics.

Practical recommendations

For particle MCMC model comparison:

Start with posterior predictive checks - do the models make qualitatively different predictions? This check applies whatever the model, and it is what the section above does.
Use an information criterion the model can supply - WAIC or LOO when each observation has a likelihood term of its own, as in the deterministic sessions, and DIC or a one-step-ahead score when a day can only be scored given the days before it, as here (see Going further).
Be cautious with Bayes factors - they are sensitive to prior choices and computationally expensive

When criteria disagree, consider what question you’re actually asking: best predictions (information criteria) or most probable model (Bayes factors)?

Learning points

Particle filters estimate likelihoods for stochastic models by simulating many trajectories (particles) and weighting them by how well they match the data
Number of particles should be calibrated so the standard deviation of the log-likelihood estimate is between 1 and 3. The required number scales with problem complexity.
Stochastic vs deterministic: stochastic models capture demographic stochasticity (important for small populations) but are much slower to fit
The model’s structure decides the criterion. In a stochastic model a day can only be scored given the days before it, which rules out WAIC and LOO, so this session uses DIC and the posterior predictive check
DIC adds twice the effective number of parameters $p_D$ to the fit at the posterior mean. SEITL and SEIT4L have the same $p_D$ , so the comparison between them comes down to fit
Posterior predictive checks apply whatever the model is, and the latent path has to come from the model for the check to have anything to fail at

Going further

Multiple chains and diagnostics

The saved chains are four chains run side by side, which is where their R̂ comes from. Turing will run several for you in one call:

# Run multiple chains in parallel
chain_multi = sample(
    model_pmmh,
    externalsampler(AdvancedMH.RobustAdaptiveMetropolis()),
    MCMCThreads(),
    10000,
    4;
    num_warmup=10000,
    check_model=false
)

# Check R̂ and ESS statistics
summarystats(chain_multi)

Use the same diagnostics (R̂ < 1.1, ESS > 100) as for standard MCMC. A single chain still reports an R̂, by splitting itself in half, so seeing one is no guarantee that several chains were run. One chain can look perfectly healthy on its own trace and still have missed a region of the posterior, which is the failure only several chains started apart can show you.

Alternative: Manual proposal tuning

If you need more control, you can use MH() with a covariance matrix instead of RAM. This is the same approach we used in the MCMC session:

# Manual covariance matrix (6 parameters, so 6x6 matrix)
# Diagonal entries are proposal variances, off-diagonal are covariances
Σ = [4.0  0.0  0.0  0.0  0.0  0.0;   # R_0
     0.0  0.04 0.0  0.0  0.0  0.0;   # D_lat
     0.0  0.0  0.25 0.0  0.0  0.0;   # D_inf
     0.0  0.0  0.0  0.01 0.0  0.0;   # α
     0.0  0.0  0.0  0.0  4.0  0.0;   # D_imm
     0.0  0.0  0.0  0.0  0.0  0.01]  # ρ

chain = sample(model_pmmh, MH(Σ), 5000)

This requires choosing proposal variances for each parameter - too large gives low acceptance, too small gives slow exploration. In practice, RAM is easier because it learns the covariance structure automatically.

Try it and compare the acceptance rate against the one RAM reached above. Matching that by hand takes several attempts, and the six variances above are only a starting point rather than a tuned answer.

Particle Gibbs

Particle Gibbs (PG) is another PMCMC method, available in Turing via PG(n_particles).

The methodological difference: PMMH uses the particle filter to estimate the marginal likelihood, then does standard Metropolis-Hastings on the parameters. PG conditions on a “reference trajectory” - one particle path from the previous iteration is retained, and the other particles are sampled conditional on it.

PG can exhibit high autocorrelation in some settings - Lindsten et al. (2014) discuss this and propose the Particle Gibbs with Ancestor Sampling (PGAS) algorithm as an improvement.

Filtered trajectories

PMCMC methods can provide filtered trajectories - stochastic trajectories from the model conditioned on both the estimated parameters and the observed data. These are useful for visualising what the underlying epidemic dynamics might have looked like, given what we observed.

The model comparison above plots these against the re-simulated replicates, using filtered_incidence from the course package. It samples one particle at the final time point, weighted by the particle weights, and reconstructs its path back through time, which gives one draw from the posterior over trajectories. Repeating across MCMC iterations gives a collection of plausible epidemic histories.

One caveat that plot does not show: within a single filter run the surviving particles share ever fewer distinct ancestors as you trace them back, so the early part of any one path is drawn from a handful of histories rather than from the full smoothing distribution. Drawing many paths from one run would therefore understate the uncertainty badly. The plot above avoids that by re-running the filter for every replicate, at a fresh posterior draw, so the paths are independent of each other.

See Doucet & Johansen (2009) for details on trajectory sampling from particle filters.

Calibrating the particle count and the proposal together

This session calibrated the particle count in one place and let RAM adapt the proposal in another, as though the two were separate problems. They are not. A rejection costs a particle filter run once the likelihood is estimated, so the count and the proposal are one calibration. Sherlock, Thiery, Roberts & Rosenthal (2015) work out the efficient combination for pseudo-marginal random walk Metropolis: a log-likelihood standard deviation near 1.8, with an acceptance rate near 7%. Both figures are limits for a large number of parameters and assume a particular form for the noise, as the 23.4% in the MCMC session does, so they are a guide rather than a setting to copy.

Hitting either target still takes a pilot run. Abaoubida, Bédard & Maire (2025) instead adapt the particle count as the chain runs, and prove that the sampler stays ergodic as long as the adaptation fades out.

SMC²

SMC² (SMC squared) uses Sequential Monte Carlo for both parameters and state trajectories, avoiding MCMC entirely. This can be more efficient for complex models but requires more careful tuning. Turing.jl supports this through the SMC sampler. Adapting the number of state particles is further along here than it is for PMMH: Botha, Kohn, South & Drovandi (2023) change it automatically as the algorithm runs, using the expected squared jumping distance of the rejuvenation moves to decide when. See Chopin, Jacob & Papaspiliopoulos (2013), which introduced the algorithm, for how the two layers of particles fit together.

Amortised inference

Every chain in this session answered a question about a single dataset. Fit the same model to a second island, or to a different reporting schedule, and the whole cost is paid again.

Neural posterior estimation is an advanced session that trains a neural network on simulations from this same SEIT4L model, so that afterwards a posterior costs milliseconds for any dataset the model can produce. It is approximate where PMMH is exact, and it pays off when you fit the same model repeatedly.

Other ways to compare the models

DIC is one of several criteria, and the others need either a different kind of model or more computation than this session does.

WAIC and LOO are the usual recommendation where each observation has a likelihood term of its own, as in the deterministic sessions. ParetoSmooth.jl computes LOO.
A one-step-ahead score judges each day against a prediction made from the days before it (Millar 2015), which suits a time series. The particle filter builds its likelihood estimate from exactly these one-step-ahead predictions, so it computes them along the way.
The marginal likelihood, for a Bayes factor. The particle filter’s estimate of $p(y \mid \theta)$ is unbiased, so averaging it over draws from the prior gives an unbiased, if very noisy, estimate of the marginal likelihood, and SMC² gives one more efficiently. Either way the answer depends on the priors.

Next session

In Observation Models and Distance Metrics, we’ll step back and examine a fundamental question: why do we choose Poisson likelihoods, and what are the alternatives?

References

Andrieu et al. (2010). Particle Markov chain Monte Carlo methods. JRSS B 72(3), 269-342. (The original PMCMC paper)
Pitt et al. (2012). On some properties of Markov chain Monte Carlo simulation methods based on the particle filter. Journal of Econometrics 171(2), 134-151. (Particle calibration: target log-likelihood SD of 1-3)
Vihola (2012). Robust adaptive Metropolis algorithm with coerced acceptance rate. Statistics and Computing 22, 997-1008. (RAM algorithm used in this session)
Roberts & Rosenthal (2009). Examples of Adaptive MCMC. J. Computational and Graphical Statistics 18(2), 349-367. (Adaptive MCMC theory)
Doucet & Johansen (2009). A tutorial on particle filtering and smoothing. Handbook of Nonlinear Filtering 12, 656-704. (Trajectory sampling)
Doucet et al. (2015). Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator. Biometrika 102(2), 295-313.
Doucet et al. (2001). Sequential Monte Carlo Methods in Practice. Springer. (Comprehensive textbook)
Arulampalam et al. (2002). A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal Processing 50(2), 174-188.
Turing.jl documentation - Probabilistic programming in Julia

References

Celeux, Gilles, Florence Forbes, Christian P. Robert, and D. Mike Titterington. 2006. “Deviance Information Criteria for Missing Data Models.” Bayesian Analysis 1 (4): 651–73. https://doi.org/10.1214/06-BA122.

Millar, Russell B. 2015. “A One-Step-Ahead Pseudo-DIC for Comparison of Bayesian State-Space Models.” Biometrics 71 (1): 201–6. https://doi.org/10.1111/biom.12237.