Sequential Monte Carlo and particle filters
We know how to get a posterior when the model is deterministic:
\[p(\theta \mid y) \propto p(y \mid \theta)\, p(\theta)\]
One \(\theta\) gives one trajectory, so the likelihood is a single calculation: solve the ODE, evaluate the observation model, done.
Stochastic models break this.
The model has a latent state \(x\) — the trajectory it actually took — which we do not observe. To get the likelihood of \(\theta\) we must marginalise over every trajectory it could have taken:
\[p(y \mid \theta) = \sum_{x} p(y \mid x, \theta)\, p(x \mid \theta)\]
For a deterministic model, \(\theta\) determines the trajectory exactly:
\[p(y \mid \theta) = \sum_{x} p(y \mid x, \theta) \times \mathbb{1}_{x = f(\theta)} = p\big(y \mid x = f(\theta), \theta\big)\]
The sum collapses to a single term. That is why everything has worked so far.
\[p(y \mid \theta) = \sum_{x} p(y \mid x, \theta)\, p(x \mid \theta)\]
Now \(x\) is no longer known, and the number of possible trajectories is astronomical — every combination of every event at every time.
We cannot enumerate them.
Approximate the sum by a sample of \(J\) particles, each one a simulated trajectory:
\[p(y \mid \theta) \approx \frac{1}{J}\sum_{j=1}^{J} p(y \mid x_j, \theta)\]
This is sequential Monte Carlo, better known as particle filtering.
Run the particles forward through time, one observation at a time:
Initialise — draw \(J\) particles from the initial state distribution, each with weight \(1/J\)
Propagate — simulate each particle forward to the next observation time
Weight — score each particle by how well it explains the observation, \(w_j = p(y_t \mid x_j, \theta)\)
Resample — draw a new set of particles in proportion to their weights, so good ones are copied and bad ones die
Repeat to the end of the data. The average weight along the way gives an estimate of \(p(y \mid \theta)\).
Without resampling, after a few time steps almost all the weight sits on one or two particles and the rest contribute nothing.
This is particle degeneracy, and it is the central practical problem. Resampling concentrates effort where the data say the trajectory actually went.
The cost is depletion: repeatedly copying survivors means the particles share ancestry, so diversity falls over time.
The likelihood estimate is unbiased but stochastic: run it twice with the same \(\theta\) and you get two different numbers. That has consequences we deal with in the next session.
This is the actual implementation you will use in the practical — propagate, weight, resample, accumulate:
function particle_filter_seit4l(θ, obs, n_particles::Integer;
init_state=[279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0])
n_particles > 0 || throw(ArgumentError("n_particles must be > 0"))
n_obs = length(obs)
ρ = θ[:ρ]
particles = [copy(init_state) for _ in 1:n_particles]
log_lik = 0.0
for t in 1:n_obs
# Propagate particles
inc = [gillespie_step_seit4l!(particles[i], θ) for i in 1:n_particles]
# Weight by observation likelihood
log_w = [logpdf(Poisson(max(ρ * inc[i], 1e-10)), obs[t]) for i in 1:n_particles]
# Log-sum-exp trick for numerical stability
max_lw = maximum(log_w)
w = exp.(log_w .- max_lw)
log_lik += max_lw + log(mean(w))
# Normalize weights
w ./= sum(w)
# Resample if ESS too low
ess = 1.0 / sum(w.^2)
if ess < n_particles / 2
idx = wsample(1:n_particles, Weights(w), n_particles)
particles = [copy(particles[i]) for i in idx]
end
end
return log_lik
endsrc/seit4l_bootstrap_filter.jl
In the practical you will
Fitting stochastic models