Particle filters

Estimated time: 105 minutes (includes computation time for calibration)

Introduction

Where we are: You’ve completed the core workflow (Introduction → MCMC → Diagnostics → Model Checking) and seen the SEITL model. Now we tackle stochastic models, where the same parameters produce different trajectories each run.

In the previous sessions, we fitted deterministic models using MCMC. For deterministic models, computing the likelihood is straightforward: given parameters, we solve the ODEs to get a trajectory, then evaluate the probability of the data given that trajectory.

But what about stochastic models? In the Structure and stochasticity session, we saw that the same parameters can produce very different trajectories due to demographic stochasticity. There’s no single “expected” trajectory to compare against the data.

This session introduces particle filters — the computational technique that makes Bayesian inference possible for stochastic models.

Objectives

In this session you will:

Understand why stochastic models need different inference methods
Learn the core particle filter algorithm: propagate, weight, resample
Implement the SEIT4L model using the state-space model interface
Calibrate the number of particles needed for reliable inference

Reminder: Stochastic variability

Particle filters use stochastic simulation, so your results will differ slightly from others’ - and from repeated runs. The patterns (e.g., how ESS changes with particle count) should be similar, but exact numbers will vary. This is expected.

Computational requirements

Particle filters are computationally intensive. For this session:

Memory: ~2-4 GB RAM for 256 particles (scales linearly with particle count)
Runtime: Single likelihood evaluation takes 0.1-0.5 seconds with 100 particles
Calibration section: Running 30 replicates at each of eight particle counts takes 1-2 minutes

If you’re on a slower machine, reduce particle counts in the calibration exercises (e.g., use 50, 100, 200 instead of 64, 128, 256, 512).

Setup

Source file

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

Load libraries

using Random ## for random numbers
using Distributions ## for probability distributions
using DataFrames ## for data frames
using Plots ## for plots
using StatsPlots ## for statistical plots
using CSV ## for CSV file reading
using StatsBase ## for statistical functions
using SSMProblems ## for state-space model interface
using GeneralisedFilters ## for particle filtering
using DrWatson

Load the course package

This loads the MFIIDD course package, which provides the helper functions used throughout — including the SEIT4L jump simulator and the SSMProblems interface, both shown in the main text below.

using MFIIDD

Initialisation

We set a random seed for reproducibility.

Random.seed!(1234)

TaskLocalRNG()

In your own work, you should run prior and posterior predictive checks when fitting models. We skip them here to focus on the particle filter mechanics.

Why deterministic likelihood fails for stochastic models

First, consider the problem. With a deterministic model, the trajectory is fully determined by the parameters:

$\text{Parameters} \rightarrow \text{Unique trajectory} \rightarrow P(\text{data}|\text{trajectory})$

With a stochastic model, the same parameters can produce infinitely many different trajectories:

$\text{Parameters} \rightarrow \text{Many possible trajectories} \rightarrow \text{???}$

The likelihood we need is $P(\text{data}|\theta)$ , which marginalises over all possible trajectories:

$P(\text{data}|\theta) = \int P(\text{data}|\text{trajectory}) \times P(\text{trajectory}|\theta) \, d(\text{trajectory})$

This integral is intractable - we can’t enumerate all possible trajectories. Particle filters give us a Monte Carlo approximation.

Load data and define SEIT4L dynamics

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

# Parameters (reasonable values from literature)
θ_example = Dict(
    :R_0 => 7.0,
    :D_lat => 1.0,
    :D_inf => 4.0,
    :α => 0.5,
    :D_imm => 10.0,
    :ρ => 0.65
)

# 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

The jump simulator

Why SEIT4L instead of SEITL?

From here on, we use the SEIT4L model (Erlang-4 temporary immunity) rather than the basic SEITL (exponential temporary immunity). As we saw in the Structure and stochasticity session, the Erlang distribution gives more realistic timing for immunity loss. The original Camacho et al. (2011) paper found that SEIT4L fits the Tristan da Cunha data better. Everything below works identically with SEITL — just with fewer T sub-compartments.

The particle filter needs to simulate many stochastic trajectories. The simulation is the one introduced in the Structure and stochasticity session. The model is a list of transitions and their rates, which JumpProcesses.jl samples exactly. The transitions below come from the package source, so the page and the code cannot drift apart. The number of temporary immunity stages is read from the state, so this one description serves SEITL and SEIT4L alike:

The transitions (src/seitl_jumps.jl)

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

Key aspects:

Stoichiometry: which compartments each transition takes from and gives to
Rates: computed from the current state and the parameters
Incidence: Inc rides in the state as a slot no rate reads, so the E→I transition increments it and daily incidence is read off like a compartment
Any number of stages: the chain of T compartments is built from k, so nothing here is specific to four

How particle filters work

The particle filter approximates the likelihood integral by:

Propagate: Simulate many independent trajectories (“particles”) forward through time
Weight: At each observation time, assign weights based on how well each particle matches the data
Resample: Replace low-weight particles with copies of high-weight particles

We can visualise this with our SEIT4L model:

# Propagate 20 particles for 5 days
n_particles = 20
particles = [Float64[init_state[:S], init_state[:E], init_state[:I],
                     init_state[:T1], init_state[:T2], init_state[:T3],
                     init_state[:T4], init_state[:L]] for _ in 1:n_particles]

# Store the daily incidence for days 1-5 of each particle
trajectories = [zeros(5) for _ in 1:n_particles]

# Simulate 5 days
## Built once and reused: the transitions and rates do not change between steps
problem = seitl_jump_problem(θ_example, particles[1])

for day in 1:5
    for i in 1:n_particles
        ## advances the particle in place and returns the day's incidence
        trajectories[i][day] = seitl_jump_step!(problem, particles[i], 1.0)
    end
end

# For weighting step later. The particles were stepped five times from the initial
# state, so `trajectories[i][k]` is each particle's incidence on day k, and the
# observation to weight against is the one recorded on day 5.
day_5_predictions = [trajectories[i][5] for i in 1:n_particles]
day_5_observed = flu_tdc.obs[5]

# Plot
p = plot(title="Particle Filter: Propagation Step",
         xlabel="Time (days)", ylabel="Daily incidence")

# Particle trajectories
for i in 1:n_particles
    plot!(p, 1:5, trajectories[i], alpha=0.4, color=:steelblue, linewidth=1, label="")
end

# Observations up to day 4 in grey, and the day 5 observation the particles are
# about to be weighted against as an open circle
scatter!(p, flu_tdc.time[1:4], flu_tdc.obs[1:4],
         label="Earlier observations", markersize=4, color=:gray)
scatter!(p, [flu_tdc.time[5]], [flu_tdc.obs[5]],
         label="Current target (day 5)", markersize=7, markershape=:circle,
         markerstrokecolor=:red, markerstrokewidth=2, color=:white)

plot!(p, legend=:topright)

Each blue line is one particle trajectory. The particles spread out because each follows different random events. The open circle is the observation recorded on day 5, and it is the one used to weight the particles below.

The weighting step

Particles that predict incidence close to the observed value get higher weights:

# Compute weights (Poisson likelihood)
log_w = [logpdf(Poisson(max(θ_example[:ρ] * pred, 0.01)), day_5_observed)
         for pred in day_5_predictions]

# Normalise
particle_weights = exp.(log_w .- maximum(log_w))
particle_weights ./= sum(particle_weights)

# Visualise
p = bar(1:n_particles, particle_weights,
        xlabel="Particle", ylabel="Weight",
        title="Particle Weights at Day 5 (obs = $day_5_observed)",
        legend=false, alpha=0.7)

# Add predictions as labels
for (i, pred) in enumerate(day_5_predictions)
    annotate!(p, i, particle_weights[i] + 0.02, text("$(round(Int, pred))", 7))
end

plot!(p)

Particles predicting values closer to the observation get higher weights.

Code explainer

exp.(log_w .- maximum(log_w)) is the log-sum-exp trick for numerical stability. The raw log weights might be very negative (e.g. -500), and exp(-500) underflows to zero. Subtracting the maximum first shifts the largest value to 0, so at least one weight is exp(0) = 1, and the others are smaller but non-zero. Since we normalise by the sum afterwards, the shift cancels out — the relative weights are unchanged.
particle_weights ./= sum(particle_weights) is shorthand for particle_weights .= particle_weights ./ sum(particle_weights) — normalise in place so the weights sum to 1.

The resampling step

When weights become too uneven, we resample - drawing particles with replacement according to their weights:

# Resample particles
n_effective = 1.0 / sum(particle_weights.^2)  # Effective sample size
println("Effective sample size: $(round(n_effective, digits=1)) out of $n_particles")

# Draw new particles
resampled_idx = sample(1:n_particles, Weights(particle_weights), n_particles)

# Count how many times each original particle was selected
counts = [count(==(i), resampled_idx) for i in 1:n_particles]

p = bar(1:n_particles, counts,
        xlabel="Original particle", ylabel="Copies after resampling",
        title="Resampling: High-weight particles duplicated",
        legend=false, alpha=0.7, color=:green)
plot!(p)

Effective sample size: 5.5 out of 20

High-weight particles are duplicated, low-weight particles are discarded. This focuses computational effort on plausible trajectories.

Code explainer

1.0 / sum(particle_weights.^2) is the effective sample size (ESS). When all weights are equal, ESS equals the number of particles; when one particle dominates, ESS approaches 1. It tells you how many particles are effectively contributing to the estimate — a common trigger for deciding when to resample.
sample(1:n_particles, Weights(particle_weights), n_particles) draws n_particles indices from 1:n_particles with replacement, where each index is chosen with probability proportional to its weight. Weights(...) is a wrapper from StatsBase.jl that tells sample to interpret the vector as sampling weights.
count(==(i), resampled_idx) counts how many times each original particle was drawn. ==(i) is partial application: it creates a function “does x equal i?” that count then applies to every element of resampled_idx. This is equivalent to sum(resampled_idx .== i) but reads more naturally.

Exercise: Understand the particle filter steps

Look at the three plots above (propagation, weighting, resampling).

In the propagation plot, why do the particle trajectories spread apart over time?
In the weighting plot, particles with predicted incidence closer to the observed value get higher weights. Why does Poisson likelihood give this behaviour?
After resampling, some particles are duplicated many times while others disappear entirely. Why is this beneficial for estimating the likelihood?

Solution

Why they spread. Each particle takes its own random draws at each step. Two particles starting from the same state see different infection and recovery events, so they diverge, and the divergence compounds. That spread is not noise in the algorithm — it is the model’s genuine uncertainty about what the epidemic did, which is exactly the thing a deterministic model cannot represent.

Why Poisson weights behave that way. The weight is $p(y_t \mid x_t)$ , the Poisson probability of the observed count given the particle’s predicted incidence. Poisson probability falls off fast as the observation moves away from the mean, so a particle predicting 12 when 11 were observed scores far higher than one predicting 40. The weights are the likelihood of each particle’s guess, no more than that.

Why resampling helps. Without it, after a few steps almost all the weight sits on one or two particles and the rest contribute nothing while still costing computation — degeneracy. Resampling spends your fixed budget of particles on the regions of state space that are actually consistent with the data so far, so the next propagation step explores where it matters. You are not throwing information away: the discarded particles had negligible weight, and their contribution is already recorded in the running likelihood estimate.

Why resampling is essential: particle degeneracy

Without resampling, particle filters suffer from degeneracy: after several time steps, almost all weight concentrates on a single particle. This happens because multiplicative weights compound - a particle that’s slightly better at time 1 gets a slightly higher weight, making it even more likely to dominate at time 2, and so on.

Imagine 100 particles after 10 time steps without resampling: - 1 particle has weight 0.99 - 99 particles share weight 0.01

This means we’re effectively using just 1 trajectory to estimate the likelihood - a very poor approximation. Resampling fixes this by redistributing particles to maintain diversity, keeping the effective sample size high.

Accumulating the likelihood

The log-likelihood is estimated by summing the log of mean weights at each time step:

$\log \hat{P}(\text{data}|\theta) = \sum_{t=1}^{T} \log \left( \frac{1}{N} \sum_{i=1}^{N} w_t^{(i)} \right)$

This is unbiased in expectation, meaning if we run the particle filter many times, the average will converge to the true likelihood.

Implementing the state-space model interface

For production use, we define the model using the SSMProblems.jl interface. This separates the model from the filtering algorithm, so you can swap in a different filter without rewriting the model, within the limits of what the model can do: Going further below sets out which of the family a simulate-only model can actually use. It also integrates with Turing.jl for the Particle MCMC in the next session.

A state-space model has three components:

Initial state distribution (SEITLInitial) - where does the system start?
Transition dynamics (SEITLDynamics) - how does the state evolve? (samples the transitions above)
Observation process (PoissonObservation) - how do observations relate to the state?

Here is the full interface implementation (defined in src/seitl_ssm_interface.jl). Each type carries a docstring spelling out its state vector and the method it implements — useful for reference once you’ve read the prose above.

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

Code explainer

A few Julia features worth unpacking:

struct SEITLDynamics <: SSMProblems.LatentDynamics defines a new type — a small container that holds some data (here, the parameter dictionary θ). The <: means “is a subtype of”. SSMProblems.LatentDynamics is an abstract type: it has no fields of its own, and exists purely to mark “this thing plays the role of a latent-dynamics model”. By declaring our struct a subtype, we tell SSMProblems.jl that our type fills that role. The same pattern applies to SEITLInitial <: SSMProblems.StatePrior and PoissonObservation <: SSMProblems.ObservationProcess.
Inside the particle filter, SSMProblems.jl calls simulate(rng, dyn, ...) and distribution(obs, ...) on whatever types it’s given. The dyn::SEITLDynamics annotation on our function definition tells Julia “use this version when dyn is one of ours”, so our code gets plugged in automatically. We implement simulate for the latent dynamics (how the state evolves one step) and for the initial state (how to draw the first state). For the observation model we instead implement distribution, which returns a Distributions.jl distribution describing how observations arise from a given state.
rng::AbstractRNG is Julia’s abstract random number generator type. Any concrete RNG (the default one, a seeded Mersenne twister, etc.) counts. Threading the rng through every simulation step is what lets the particle filter control randomness — set the seed once and you get reproducible particles.
kwargs... at the end of a function signature means “accept any additional keyword arguments”. We don’t use them here, but SSMProblems.jl may pass extra keyword arguments depending on the filter being used, so our methods need to accept them without complaining.
vcat(prior.init_state, 0.0) appends 0.0 to the compartments to make the state one element longer. That last slot tracks daily incidence, because the observation process needs to see it. collect(prev_state[1:(end - 1)]) takes everything except that slot and copies it into a plain Vector — some inputs arrive as views or other array-like types, and the simulator wants a concrete vector. Writing it as end - 1 lets the same code run SEITL and SEIT4L.

How the pieces fit together:

The state vector is [S, E, I, T_1 … T_k, L, daily_inc] - the last element tracks daily incidence for the observation process, and SEIT4L is the k = 4 case
SEITLDynamics advances the state one day by sampling the transitions
PoissonObservation assumes observed cases follow Poisson(ρ × true_incidence)

Running the bootstrap particle filter

The run_particle_filter function assembles the SSM components and runs the bootstrap filter:

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

This ties everything together: the initial state, jump dynamics, and Poisson observations form a StateSpaceModel, which the bootstrap filter (BF) processes to estimate the log-likelihood.

Code explainer

StateSpaceModel(initial, dynamics, observation) is a constructor from SSMProblems.jl that bundles the three SSM components into a single object. Constructors in Julia look like function calls with the type’s name — there’s no explicit new keyword as in Python or Java.
BF(n_particles) constructs a Bootstrap Filter object from GeneralisedFilters.jl. The bootstrap filter is the simplest particle filter variant: it propagates particles using the model’s own dynamics, weights them by the observation likelihood, and resamples. A different filter can be chosen by changing this argument, although not all models can use all filters. This is further explored in the Going further section below.
GeneralisedFilters.filter(rng, model, algo, obs) returns a tuple of (particles_at_final_time, log_likelihood). We ignore the particles with _ and keep only the log-likelihood estimate, which is what MCMC will eventually need.
Random.default_rng() returns Julia’s shared global RNG. If you seed it with Random.seed!(n) at the start of the session, every call that uses this RNG becomes reproducible.

Testing it:

# Run with 100 particles
ll = run_particle_filter(θ_example, flu_tdc.obs, 100)
println("Log-likelihood estimate: $(round(ll, digits=1))")

Log-likelihood estimate: -120.2

Calibrating particle count

A key question: how many particles do we need?

Too few particles produce noisy likelihood estimates, causing poor MCMC mixing
Too many particles waste computation

The principled approach targets a specific variance of the log-likelihood estimate. Pitt et al. (2012) and Doucet et al. (2015) showed that MCMC efficiency is maximised when the standard deviation of the log-likelihood estimate is between 1 and 3.

Why is some noise acceptable? MCMC with a noisy likelihood (called “pseudo-marginal MCMC”) still targets the correct posterior — the noise slows down mixing but doesn’t introduce bias. An SD of 1-3 is a sweet spot: low enough for reasonable MCMC acceptance rates, high enough to avoid wasting computation on unnecessary precision.

Calibrating (this takes 1-2 minutes to run all particle counts):

# Reference parameters (near what we expect the posterior mode to be)
θ_cal = Dict(:R_0 => 6.0, :D_lat => 1.3, :D_inf => 2.0,
             :α => 0.5, :D_imm => 10.5, :ρ => 0.7)

# Test different particle counts
particle_counts = 2 .^ (2:9)  # 4, 8, 16, ..., 512
n_reps = 200  ## the spread is itself an estimate, and 30 draws is too few to pin it down

calibration = DataFrame(n_particles=Int[], mean_ll=Float64[],
                        sd_ll=Float64[], time_per_run=Float64[])

for n in particle_counts
    t0 = time()
    lls = Vector{Float64}(undef, n_reps)
    Threads.@threads for i in 1:n_reps
        lls[i] = run_particle_filter(θ_cal, flu_tdc.obs, n)
    end
    elapsed = (time() - t0) / n_reps

    push!(calibration, (n, mean(lls), std(lls), elapsed))
end

calibration

8×4 DataFrame

Row	n_particles	mean_ll	sd_ll	time_per_run
	Int64	Float64	Float64	Float64
1	4	-187.585	522.101	0.00200807
2	8	-209.409	758.151	0.00206267
3	16	-124.398	9.47103	0.00334324
4	32	-120.26	7.6607	0.00644618
5	64	-116.855	3.42506	0.0128111
6	128	-116.116	0.842228	0.0256849
7	256	-115.998	0.577563	0.0521487
8	512	-115.902	0.36648	0.108815

Code explainer

Threads.@threads — runs loop iterations in parallel across available threads. The replicates are independent, so this is safe and gives a roughly linear speedup.
Vector{Float64}(undef, n_reps) — pre-allocates the results vector so each thread writes to its own index (lls[i] = ...). We can’t use push! here because threads appending to a shared vector is a race condition.
To actually use more than one thread, start Julia with julia --threads=auto (or set JULIA_NUM_THREADS). Without that, @threads silently runs serially.

Plotting code

# Plot calibration results
p1 = plot(calibration.n_particles, calibration.mean_ll,
    xlabel="Number of particles", ylabel="Mean log-likelihood",
    title="Bias", legend=false, marker=:circle, linewidth=2, xscale=:log2)

p2 = plot(calibration.n_particles, calibration.sd_ll,
    xlabel="Number of particles", ylabel="SD of log-likelihood",
    title="Variance (target: 1-3)", legend=false, marker=:circle,
    linewidth=2, xscale=:log2)
hline!(p2, [1, 3], color=:grey, linestyle=:dot, alpha=0.7)

p3 = plot(calibration.n_particles, calibration.time_per_run,
    xlabel="Number of particles", ylabel="Time (seconds)",
    title="Computation time", legend=false, marker=:circle,
    linewidth=2, xscale=:log2)

plot(p1, p2, p3, layout=(1, 3), size=(900, 300))

Look at the variance plot (middle). We want the SD to be in the range 1-3. Read off the smallest particle count that achieves this.

Typical results

In our runs, we typically see:

4-16 particles: SD in the tens or worse, far too noisy for reliable inference
32 particles: SD around 7, still outside the target range
64 particles: SD around 3, at the edge of the target range
128-512 particles: SD around 1 or below, falling slowly while the cost keeps rising

Your numbers will differ, and at the low counts they differ a lot, which is why the cell above averages over 200 runs rather than a handful.

Exercise: Choose the number of particles

Based on the calibration plots above:

What particle count achieves SD of log-likelihood between 1 and 3?
What’s the trade-off between 128 and 256 particles for this problem?
Why might you want to use more particles than the minimum needed?

What to look for

Read the first answer off your own calibration plot rather than taking a number from here — it depends on the parameters you calibrated at, and the point of the exercise is that this is something you measure rather than guess.

128 against 256. Doubling the particles roughly doubles the cost per likelihood evaluation and cuts the standard deviation of the estimate by about $\sqrt{2}$ . In PMCMC that noise is what makes chains sticky, so you are trading time per iteration against the number of iterations you need. The usual advice is to aim for a standard deviation somewhere between 1 and 3 on the log scale, which is where the total cost tends to be lowest.

Why go above the minimum. The noise depends on where in parameter space you are, and a count calibrated at one parameter set will be too few elsewhere. Elsewhere is precisely where a chain gets stuck and stays stuck, so taking some headroom is cheaper than diagnosing a chain that will not move.

Exercise: Sensitivity to parameters

The particle filter likelihood estimate depends on the parameters. Explore how the log-likelihood changes as you vary R₀:

# Vary R_0 and compute log-likelihoods
R0_values = 3.0:1.0:12.0
n_particles = 128  # Use your calibrated value

ll_results = Float64[]
for R0 in R0_values
    θ_test = Dict(:R_0 => R0, :D_lat => 1.3, :D_inf => 2.0,
                  :α => 0.5, :D_imm => 10.5, :ρ => 0.7)
    ll = run_particle_filter(θ_test, flu_tdc.obs, n_particles)
    push!(ll_results, ll)
end

plot(R0_values, ll_results, xlabel="R_0", ylabel="Log-likelihood",
     title="Likelihood surface", marker=:circle, legend=false)

Questions:

Which value of R₀ gives the highest log-likelihood?
Run the same R₀ value multiple times. How much does the log-likelihood vary? (This is the noise from particle filtering)
How does the shape of this likelihood surface compare to what you’d expect from deterministic model fitting?

What to look for

Two things, and the second is the one that matters.

The log-likelihood should peak somewhere in the middle of the range you sweep and fall away on both sides, giving the usual shape you would get from a deterministic likelihood. That is reassuring: a noisy estimate still carries the signal.

But run the same $R_0$ twice and you will not get the same number. The curve is noisy, and near the peak — where it is flattest — the noise can be comparable to the difference between neighbouring parameter values. That is the whole difficulty of PMCMC in one plot: the sampler has to decide between two parameters using numbers that wobble by a similar amount as the difference it is trying to detect. More particles shrink the wobble, at a cost per evaluation.

If you want to see it plainly, compute the likelihood at a single $R_0$ ten times over and look at the spread.

Learning points

Stochastic models require particle filters because the likelihood marginalises over all possible trajectories
Particle filters work by: propagate particles forward, weight by observation likelihood, resample when weights become uneven
SSMProblems.jl interface separates model specification from filtering algorithms
Calibrate particle count to achieve log-likelihood SD between 1 and 3 (too few = noisy, too many = slow)
The likelihood estimate is unbiased but noisy - this noise slows MCMC but doesn’t bias posteriors
Next step: The noisy likelihood estimate slots directly into the Metropolis-Hastings acceptance ratio — wherever we previously computed log_likelihood(θ), we now use run_particle_filter(θ, data, n_particles) instead

Going further

GeneralisedFilters.jl provides a family of algorithms behind the interface used above. Which filter can be used depends on the simulation model used. A jump simulator can be sampled from, but its transition density cannot be evaluated, and the dynamics defined above implement simulate and nothing else.

Why the bootstrap filter. BF is not a separate algorithm from the general particle filter PF. It is PF with the proposal set to the model’s own dynamics, so the proposal and the dynamics cancel out of the importance weight and the weights reduce to the observation density alone. That cancellation is also what makes it the filter to reach for here, because it never has to evaluate the transition density. A proposal that looks at the observation before propagating would spread the weights less, but forming its importance weight needs that density, so a guided proposal stops at a MethodError on SSMProblems.distribution. AuxiliaryParticleFilter reweights by a look-ahead score instead. How that score is formed determines whether it can be used: MeanPredictive and ModePredictive evaluate the density and fail the same way, while DrawPredictive draws a step forward by simulation and runs here.

What you can vary. Resampling only duplicates and discards particles that already exist, so it proposes no new state, needs no transition density, and this axis stays open. BF takes a resampler and a threshold. Resampling cures the degeneracy described above, but it is a random step in its own right and adds variance of its own. BF(n_particles) resamples at every step, which is what every result in this session and the next one uses. A threshold of 0.5 resamples only once the effective sample size falls below half the particles, which is the common alternative. At the parameters the calibration cell uses, resampling only below the threshold brings the spread at 64 particles down to roughly what doubling the particle count gives, and at about half the cost of doubling. Which setting wins depends on the model and the particle count, so it is worth measuring on your own problem rather than carrying either default across. Substituting BF(n_particles; threshold = 0.5, resampler = GeneralisedFilters.Systematic()) for BF(n_particles) in the runner turns that argument into something you can measure with the calibration cell, qualified because the runner imports BF and not the resamplers. Systematic is already the default, so that substitution moves the threshold and nothing else. The resampler is a second experiment and not a free one: against the default, Multinomial roughly doubles the spread, so vary one at a time or the two effects will be read as one. threshold = 0.0, which never resamples at all, puts the degeneracy back so you can see what it costs.

Exploiting structure. KalmanFilter computes the likelihood exactly for linear Gaussian models, returning the same number every run rather than a noisy estimate. RBPF pairs an outer particle filter with an exact inner filter, spending particles only on the part of a model that is not linear and Gaussian. Neither fits a jump compartmental model, but both are worth knowing before reaching for particles.

Filters inside other algorithms. A filter can be conditioned on a reference trajectory through ref_state, which is what particle Gibbs needs, and that much works with a simulate-only model. Ancestor sampling does not, because it scores the reference path with the transition density that is missing here. The next session takes the other route, feeding the likelihood estimate to a Metropolis-Hastings sampler.

The GeneralisedFilters.jl documentation lists the resamplers and the rest of the family.

References

Doucet & Johansen (2009). A tutorial on particle filtering and smoothing. Handbook of Nonlinear Filtering 12, 656-704. (Excellent introduction)
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. (Calibration: target SD 1-3)
Doucet et al. (2015). Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator. Biometrika 102(2), 295-313.
SSMProblems.jl documentation - State-space model interface
GeneralisedFilters.jl documentation - Particle filtering algorithms

Next session: In Particle MCMC, we’ll use the particle filter within an MCMC framework to fit the stochastic SEIT4L model to data.