Estimating the likelihood

Sequential Monte Carlo and particle filters

Where we are

We know how to get a posterior when the model is deterministic:

p(θ∣y)∝p(y∣θ)p(θ)p(\theta \mid y) \propto p(y \mid \theta)\, p(\theta)

One θ\theta gives one trajectory, so the likelihood is one calculation.

Stochastic models break this.

The marginal likelihood

The model has a latent state xx, the trajectory it actually took, which we do not observe.

p(y∣θ)=∫p(y∣x,θ)p(x∣θ)dxp(y \mid \theta) = \int p(y \mid x, \theta)\, p(x \mid \theta) \, \mathrm{d}x

Getting the likelihood of θ\theta means averaging over every trajectory the model could have taken.

The deterministic case

θ\theta determines the trajectory exactly:

p(y∣θ)=∫p(y∣x,θ)δ(x−f(θ))dx=p(y∣x=f(θ),θ)p(y \mid \theta) = \int p(y \mid x, \theta)\, \delta\big(x - f(\theta)\big) \, \mathrm{d}x = p\big(y \mid x = f(\theta), \theta\big)

The integral collapses to a single point. That is why everything has worked so far.

The stochastic case

p(y∣θ)=∫p(y∣x,θ)p(x∣θ)dxp(y \mid \theta) = \int p(y \mid x, \theta)\, p(x \mid \theta) \, \mathrm{d}x

The integral runs over every sequence of events the model could produce, and the times at which they happen.

Estimating the integral

One draw

Simulating the model draws x1∼p(x∣θ)x_1 \sim p(x \mid \theta). Scoring it against the data gives p(y∣x1,θ)p(y \mid x_1, \theta), an unbiased but very noisy estimate of the likelihood.

Many draws

p̂(y∣θ)=1J∑j=1Jp(y∣xj,θ)\hat{p}(y \mid \theta) = \frac{1}{J}\sum_{j=1}^{J} p(y \mid x_j, \theta)

Each trajectory is one way the epidemic could have gone. Its weight is the probability of the data if it had gone that way. A trajectory with a higher weight counts for more in the average.

One weight per trajectory

One weight has to cover all 59 observations. Two of the 200 hold 88% of the total.

The estimate gets noisier every day

Scored on day 1 alone, the effective sample size is close to all 200. By day 10 it is 6, and by day 30 about 1.

The estimate then rests on one or two draws, so it changes wildly from run to run.

The fall is close to a straight line on a log scale, so a bigger JJ cannot keep up.

The particle filter

The likelihood factorises

p(y∣θ)=p(y1∣θ)p(y2∣y1,θ)⋯p(yT∣y1:T−1,θ)p(y \mid \theta) = p(y_1 \mid \theta)\; p(y_2 \mid y_1, \theta) \cdots p(y_T \mid y_{1:T-1}, \theta)

Each factor asks how well the model predicted one day, given everything before it.

We track the possible trajectories with a set of JJ particles.

Initialise

JJ particles start from the same state, each with weight 1/J1/J.

Propagate

Every particle takes one step of the model.

Weigh

The weight is the probability of today’s count under that particle, wt,j=p(yt∣xt,j,θ)w_{t,j} = p(y_t \mid x_{t,j}, \theta).

Weighting alone changes nothing

Let the weights accumulate across days, so each particle carries the product of its own weights:

1J∑j∏twt,j=1J∑jp(y∣xj,θ)\frac{1}{J}\sum_{j} \prod_{t} w_{t,j} = \frac{1}{J}\sum_{j} p(y \mid x_j, \theta)

That is the estimator we started with. The factorisation gives us a point between days at which to discard particles.

Otherwise all 200 are simulated to day 59, including the 198 that contribute almost nothing.

Resample

Draw JJ particles from the current ones in proportion to their weights. Every weight resets to 1/J1/J.

Five days of propagate, weigh, resample

Four particles. Red diamonds are the data, yellow circles the particles, sized by weight, struck through when nobody resampled them.

The estimator

Average the weights each day, before resampling:

p̂(y∣θ)=∏t=1T1J∑j=1Jwt,j\hat{p}(y \mid \theta) = \prod_{t=1}^{T} \frac{1}{J}\sum_{j=1}^{J} w_{t,j}

The average runs over all JJ particles, including the ones resampling is about to discard. That is why resampling cannot inflate the estimate.

This is sequential Monte Carlo.

Particle degeneracy

Without resampling the effective sample size is 6 by day 10, and about 1 from day 30 on.

Depletion

Follow the arrows back from day 5. The survivors all descend from a few of the day 1 particles.

Early states end up represented by fewer and fewer distinct particles, which limits what the filter can say about the start of the epidemic.

How many particles?

  • Too few and the likelihood estimate is noisy, which will wreck the sampler that uses it.
  • Too many and every likelihood evaluation becomes expensive.

The estimate is unbiased but stochastic. Run it twice at the same θ\theta and you get two different numbers.

Measure the spread of log⁡p̂\log \hat{p} at one θ\theta and pick JJ from that. You do this in the practical.

The whole thing, in code

function bootstrap_filter(rng, θ, obs, n_particles; init_state)
    particles = [copy(init_state) for _ in 1:n_particles]
    problem = seitl_jump_problem(θ, init_state; rng = rng)
    loglik = 0.0

    for y in obs
        ## propagate: run every particle forward one day
        incidence = map(1:n_particles) do i
            seitl_jump_step!(problem, particles[i], 1.0)
        end

        ## weight: how plausible is today's observation under each particle?
        log_w = logpdf.(Poisson.(θ[:ρ] .* incidence), y)

        ## accumulate: the mean weight estimates today's likelihood
        loglik += logsumexp(log_w) - log(n_particles)

        ## resample: keep the particles that explained the data
        w = exp.(log_w .- maximum(log_w))
        particles = particles[rand(rng, Categorical(w ./ sum(w)), n_particles)]
    end

    return loglik
end

In the practical you read this one line by line, then swap it for the bootstrap filter in GeneralisedFilters.jl behind the same interface.

Your Turn

In the practical you will

  • see why a deterministic likelihood fails for a stochastic model
  • watch particles propagate, get weighted, and be resampled
  • read that filter line by line and work out what each step does
  • explore how the number of particles affects the likelihood estimate

Return to the session