Approximate Bayesian computation

Inference by simulation

What if you cannot write down a likelihood?

Imagine a model you can simulate but cannot score. Run it and out comes a dataset. Ask how probable your dataset was under it, and there is no formula to evaluate.

Reporting delays, weekend effects, behaviour that responds to the epidemic: each is simple to simulate, and together they have no density you can write.

The likelihood is still there. The simulator is p(y∣θ)p(y \mid \theta) — you simply cannot evaluate it.

Draw, simulate, compare

Draw θ\theta from the prior. Simulate a dataset. Keep θ\theta if the simulation came out exactly equal to the data.

Nothing has been approximated. The θ\theta you keep are draws from p(θ∣y)p(\theta \mid y).

Approximate Bayesian computation

For continuous data exact equality has probability zero. For our 59 counts it is possible, but so unlikely that we would wait effectively forever. Two concessions, and this is where the approximate comes from.

Close, not equal. Keep θ\theta when the distance dd between simulation and data comes in under a tolerance ϵ\epsilon.

Summaries, not the whole dataset. Over 59 days even “close” almost never happens, so compare summary statistics SS instead.

p(data∣θ)≈p(d(Sdata,Ssim(θ))<ϵ)p(\text{data} \mid \theta) \approx p\big(d(S_\text{data},\, S_{\text{sim}(\theta)}) < \epsilon\big)

Choosing summary statistics

How many. All of them must land inside ϵ\epsilon at once, so acceptance falls off fast as you add more. Three or four is usual.

Whether they move with θ\theta. One that barely responds shrinks the acceptance region and tells you nothing.

What scale they are on. Weight each by its spread under the prior predictive, or the largest one drowns out the rest.

Sufficient summary statistics

A summary is sufficient if it tells you everything the data say about θ\theta. Then p(θ∣S(y))p(\theta \mid S(y)) is p(θ∣y)p(\theta \mid y).

For a model like ours there is no such summary short of the whole dataset.

removed by
the threshold ϵ\epsilon letting ϵ→0\epsilon \to 0
the summary SS nothing

Choosing ϵ\epsilon

Too tight and you reject almost everything. The acceptance rate collapses and you wait forever for a handful of samples.

Too loose and you accept almost everything, including parameters that fit badly. In the limit you are sampling the prior.

There is no free lunch here: ϵ\epsilon trades accuracy against compute, and you cannot tell from inside the algorithm whether you have chosen well.

Rejection ABC

The simplest version, and the one you will write:

  1. Sample θ*\theta^* from the prior
  2. Simulate a dataset from the model at θ*\theta^*
  3. Compute the distance dd between simulation and data
  4. Keep θ*\theta^* if d<ϵd < \epsilon, otherwise discard it

Repeat until you have enough samples. The kept parameters approximate the posterior.

No likelihood anywhere. Just a simulator, a distance, and a threshold.

ABC-SMC

Rejection ABC draws every parameter from the prior, however hopeless. ABC-SMC walks the threshold down in stages instead.

  1. Run rejection ABC at a loose ϵ1\epsilon_1
  2. Perturb the accepted particles to propose for the next round
  3. Tighten to ϵ2<ϵ1\epsilon_2 < \epsilon_1 and repeat

Each generation starts from parameters already known to be reasonable, so a tight threshold becomes reachable.

The catch: particles need weighting to stay unbiased, which is the same propagate-weight-resample structure as the particle filter.

Your Turn

In the practical you will

  • read an ABC rejection sampler and choose its threshold
  • find a workable acceptance threshold, and see why that is the hard part
  • trace ABC-SMC and compare it against rejection sampling
  • see what a summary statistic and a distance each do to the posterior
  • compare the ABC posterior against particle MCMC on the same three parameters, and on all six

Return to the session