Particle MCMC

Fitting parameters of stochastic models

Plugging in an estimated likelihood

Where we got to

The particle filter gives us an estimate of the marginal likelihood for one value of \(\theta\):

\[\hat{p}(y \mid \theta) \approx \frac{1}{J}\sum_{j=1}^{J} p(y \mid x_j, \theta)\]

But we do not want the likelihood at one \(\theta\). We want the posterior over \(\theta\).

The obvious idea

We already have an algorithm that explores a posterior using nothing but likelihood evaluations: Metropolis-Hastings.

So: propose \(\theta^*\), run a particle filter to get \(\hat{p}(y \mid \theta^*)\), and accept or reject as usual.

\[r = \min\left(1, \frac{\hat{p}(y \mid \theta^*)\, p(\theta^*)}{\hat{p}(y \mid \theta)\, p(\theta)}\right)\]

This is particle marginal Metropolis-Hastings (PMMH).

Why should that work at all?

We are plugging a noisy, random likelihood estimate into an algorithm derived assuming the exact likelihood.

Run the filter twice at the same \(\theta\) and you get two different numbers. It looks like this should break the chain.

It does not. The chain still targets the correct posterior — exactly, not approximately.

The reason it works

The particle filter estimate is unbiased:

\[\mathbb{E}\big[\hat{p}(y \mid \theta)\big] = p(y \mid \theta)\]

Treat the random numbers used inside the filter as auxiliary variables and run Metropolis-Hastings on the enlarged space of \((\theta, \text{particles})\).

That extended chain has the right invariant distribution, and \(\theta\) marginally has the posterior we wanted.

The noise costs us efficiency, never correctness.

Andrieu et al. (2010)

The practical problem

Noise makes the chain sticky

When the filter happens to return an unusually high estimate, that value sits in the denominator of every subsequent acceptance ratio.

The chain then rejects good proposals for a long stretch, because it is comparing them against a lucky number rather than the truth.

The result is the flat plateaus you will see in the trace plots — much worse mixing than any chain we have run so far.

Choosing the number of particles

There is a genuine trade-off:

  • Too few particles: the likelihood estimate is noisy, the chain sticks, and mixing collapses.
  • Too many particles: every single MCMC iteration costs a full particle filter run, so the chain crawls.

The standard guidance is to tune \(J\) so that the standard deviation of the log likelihood estimate is roughly 1 to 3 at a representative \(\theta\).

Pitt et al. (2012); Doucet et al. (2015)

Why RAM rather than NUTS

We cannot use gradients: the particle filter involves resampling and discrete event simulation, so the likelihood is not differentiable.

That rules out NUTS and leaves us with random-walk methods, which is exactly the fallback flagged in the MCMC session.

We use robust adaptive Metropolis, so at least the proposal covariance is learned rather than hand-tuned (Vihola 2012).

What to expect when you run it

  • Acceptance rates lower than you are used to
  • Long flat stretches in the trace
  • Convergence needing far more iterations than the deterministic fits

The pre-computed chain in the practical used 500,000 samples, 50,000 burn-in, thinned by 50, and took hours. This is normal for PMCMC, and it is the price of fitting a genuinely stochastic model.

The pay-off

Everything from the deterministic part of the course now carries over to stochastic models:

  • posterior distributions over parameters
  • credible intervals
  • posterior predictive checks

And we get the latent trajectories for free: each accepted \(\theta\) comes with a sampled state path, which is state estimation alongside parameter estimation.

Comparing models, not just fitting one

Is temporary immunity better described by one compartment or four?

That is not a question about parameters. No amount of sampling SEITL tells you whether SEIT4L would have done better — you have to fit both.

Ways to compare:

  • Posterior predictive checks — which reproduces the features you care about?
  • Predictive performance on held-out data
  • The marginal likelihood, if you can estimate it — expensive, and the particle filter already gives you an estimate of it

In the practical you fit both and compare. Expect them to agree on the first wave and differ on the timing of the second.

Your Turn

In the practical you will

  • assemble the particle filter and Metropolis-Hastings into PMMH
  • run a short chain and see the poor mixing for yourself
  • vary the particle count and watch the trade-off appear
  • work with a pre-computed long chain, and compare against the deterministic fit

References

Andrieu, Christophe, Arnaud Doucet, and Roman Holenstein. 2010. “Particle Markov Chain Monte Carlo Methods.” Journal of the Royal Statistical Society Series B 72 (3): 269–342. https://doi.org/10.1111/j.1467-9868.2009.00736.x.
Doucet, Arnaud, Michael K. Pitt, George Deligiannidis, and Robert Kohn. 2015. “Efficient Implementation of Markov Chain Monte Carlo When Using an Unbiased Likelihood Estimator.” Biometrika 102 (2): 295–313. https://doi.org/10.1093/biomet/asu075.
Pitt, Michael K., Ralph dos Santos Silva, Paolo Giordani, and Robert Kohn. 2012. “On Some Properties of Markov Chain Monte Carlo Simulation Methods Based on the Particle Filter.” Journal of Econometrics 171 (2): 134–51. https://doi.org/10.1016/j.jeconom.2012.06.004.
Vihola, Matti. 2012. “Robust Adaptive Metropolis Algorithm with Coerced Acceptance Rate.” Statistics and Computing 22 (5): 997–1008. https://doi.org/10.1007/s11222-011-9269-5.

Return to the session