Approximate Bayesian computation

Inference when you cannot write down a likelihood

Estimated time: 60 minutes | Requires: Observation models

Slides

Introduction

For many models you cannot evaluate $p(y \mid \theta)$, but you can still run the model forward. Approximate Bayesian computation is built on that alone: draw parameters from the prior, simulate a dataset, and keep the parameters whose simulation came out equal to the data. Nothing is approximated at that point, and the parameters you keep are draws from the posterior.

ABC runs the simulator from beginning to end and compares the output with the data through a distance you choose. That is the whole of what it needs from a model. It never opens the simulator up, and it never asks for the probability of an observation.

Objectives

In this session you will:

1. Say what ABC needs from a model, and what it does not
2. Read an ABC rejection sampler and choose a workable acceptance threshold, which is the hard part of the method
3. Explain why ABC-MCMC stalls at a tight threshold, and how ABC-SMC reaches one
4. Compare an ABC posterior with a likelihood-based one that has the same target, and tell the approximation apart from the choice of which parameters to fix
5. Say what a summary statistic and a distance each do to the posterior, and why a tighter threshold cannot rescue a summary that is not sufficient

Setup

Source file

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

Load libraries

using Distributions
using DataFrames
using CSV
using Plots
using StatsPlots
using Turing
using Random
using LinearAlgebra
using DrWatson
using StatsBase  # For Weights in ABC-SMC
using FlexiChains ## for ess, used by the ABC-MCMC diagnostics
using MFIIDD ## the course package, which holds the SEIT4L implementation

Precompiling packages...
   1003.3 ms  ✓ QuartoNotebookWorkerTablesExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
    823.4 ms  ✓ QuartoNotebookWorkerLaTeXStringsExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   1712.8 ms  ✓ QuartoNotebookWorkerJSONExt (serial)
  1 dependency successfully precompiled in 2 seconds
Precompiling packages...
   2070.3 ms  ✓ QuartoNotebookWorkerPlotsExt (serial)
  1 dependency successfully precompiled in 2 seconds

Initialisation

We set a random seed for reproducibility.

Random.seed!(1234)

TaskLocalRNG()

Load data

We use the same influenza outbreak data from 1971:

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

# Plot
scatter(
    flu_tdc.time, flu_tdc.obs,
    xlabel = "Time (days)", ylabel = "Daily incidence",
    label = "Observed cases", markersize = 4, color = :red,
    legend = :topright
)

Daily reported cases on Tristan da Cunha in 1971. Day 1 is 13 August, when the ship arrived.

We have 59 daily observations of case counts.

We also load two particle MCMC chains, precomputed in data/, to compare the ABC posteriors against. The first is the fit from the previous session, which estimated all six SEIT4L parameters. The second estimates only the three that ABC estimates below, with the other three fixed at the same values, so it is aiming at the same posterior as ABC. It was generated by scripts/generate_pmcmc_seit4l_abc.jl.

pmcmc_chain = CSV.read(datadir("pmcmc_seit4l_chain.csv"), DataFrame) ## all six parameters
pmcmc_three = CSV.read(datadir("pmcmc_seit4l_abc_chain.csv"), DataFrame) ## as ABC
first(pmcmc_chain, 3)

3×13 DataFrame

Row	chain	R_0	D_lat	D_inf	α	D_imm	ρ	logα	η	accepted	logprior	loglikelihood	logjoint
	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Bool	Float64	Float64	Float64
1	1	6.41907	1.77069	1.60579	0.469772	10.323	0.730944	-2.50517	0.00151572	false	-8.19946	-117.275	-125.474
2	1	6.19723	1.17365	1.93196	0.422395	8.32502	0.688121	-3.68801	0.00151572	false	-8.27754	-118.715	-126.993
3	1	6.25192	1.01462	2.56485	0.371971	11.3423	0.590976	0.0	0.00151572	true	-8.11534	-118.544	-126.659

The model we fit

We keep the SEIT4L model from the particle filter and particle MCMC sessions, simulated stochastically. Run it twice at the same parameters and you get two different epidemics, so there is no single curve to hold against the data.

The particle MCMC session already fitted this model, and fitted it better than this session will. It estimated all six parameters where we estimate three, and it reached a tighter posterior from fewer simulations. On this outbreak the particle filter is the better tool, and we use the same model anyway so that you can hold ABC against an answer you already trust.

ABC earns its place where the particle filter cannot run at all. The filter weights each particle by an observation density $p(y_t \mid x_t)$, so you have to be able to write one down. It also stops the simulator at every observation, copies the state and restarts from it, so you have to be able to checkpoint the simulator. A compiled model you cannot open, someone else’s executable, or a reporting process nobody can put a density on: any of these rules the filter out and leaves ABC standing. Neither is true here, which is exactly why this is a good place to learn ABC and a poor argument for using it.

Three parameters are estimated, $R_0$, $D_\mathrm{inf}$ and $\rho$, and the rest are fixed at the values the earlier sessions used. Every free parameter has to land inside the acceptance window at the same time, so the acceptance rate falls away sharply as the number of free parameters grows.

## Fixed at the values used in the structure and stochasticity session
const D_LAT = 2.0
const ALPHA = 0.5
const D_IMM = 13.0

## The island as the ship lands, at the start of day 1: three of the five
## returning islanders are past their infectious period and start in T1, the two
## who fell ill on landing start in I. N = 284
init_test = [279.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0]
n_obs = length(flu_tdc.obs)

## One stochastic SEIT4L run, thinned to reported cases, one value per day
function simulate_seit4l_obs(R_0, D_inf, ρ)
    θ = Dict(:R_0 => R_0, :D_lat => D_LAT, :D_inf => D_inf,
             :α => ALPHA, :D_imm => D_IMM)
    state = copy(init_test)
    problem = seitl_jump_problem(θ, state)
    inc = [seitl_jump_step!(problem, state, 1.0) for _ in 1:n_obs]
    ## The floor matches the particle filter's observation model, where it keeps a
    ## particle's weight finite on a day it has no cases; here it only keeps the
    ## two models the same
    return [rand(Poisson(max(ρ * x, 1e-10))) for x in inc]
end

simulate_seit4l_obs (generic function with 1 method)

p_runs = plot(xlabel = "Time (days)", ylabel = "Cases")

for i in 1:10
    plot!(p_runs, flu_tdc.time, simulate_seit4l_obs(7.0, 2.5, 0.65),
          color = :steelblue, alpha = 0.4, label = i == 1 ? "Simulated" : "")
end

scatter!(p_runs, flu_tdc.time, flu_tdc.obs, color = :red, markersize = 4,
         label = "Observed")

Ten stochastic runs of the model at $R_0 = 7$, $D_\mathrm{inf} = 2.5$ days and $\rho = 0.65$, against the reported cases.

Ten runs, one parameter set, ten different epidemics. That spread sets the standard every distance in this session has to meet, because those ten runs are all equally plausible under those parameters. It prevents even a good parameter set from landing close to the data reliably, and the ABC-MCMC section below shows what that does to a sampler.

NoteThe likelihood fits we compare against

Throughout the session we hold the ABC posteriors against likelihood-based ones, from the two particle MCMC chains loaded with the data. The Poisson observation model is still in use, inside the filter, weighting particles one observation at a time. The ABC simulator draws its reported cases from the same Poisson distribution, so the two methods fit the same model and differ in how they use it. What has no closed form is $p(y \mid \theta)$ itself, because forming it means summing over every path the epidemic could have taken. The particle filter estimates that sum. ABC never forms it, and compares simulated output with the data through a distance.

Approximate Bayesian computation

ABC: Making the distance explicit

A parametric likelihood settles three questions for you, quietly, as a consequence of the distribution you picked. ABC puts the same three in front of you and makes you answer them yourself:

1. Summary statistics: which features of the data matter?
2. Distance: how do we measure similarity?
3. Threshold: how close is close enough?

That is the sense in which ABC makes the distance explicit. The choices were always there.

Here is a simple ABC rejection sampler from scratch:

"""
Simple ABC rejection sampler.
Samples from prior, simulates data, accepts if distance < threshold.
"""
function abc_rejection(prior_sampler, simulator, distance_fn, observed_data,
                      threshold, n_samples)
    samples = []
    n_attempts = 0
    max_attempts = n_samples * 1000  # Give up after many attempts

    while length(samples) < n_samples && n_attempts < max_attempts
        n_attempts += 1

        # Sample from prior
        θ = prior_sampler()

        # Simulate data
        simulated_data = simulator(θ)

        # Compute distance
        d = distance_fn(simulated_data, observed_data)

        # Accept if close enough
        if d < threshold
            push!(samples, θ)
        end
    end

    acceptance_rate = length(samples) / n_attempts
    println("Acceptance rate: $(round(acceptance_rate * 100, digits=2))%")
    println("Attempts needed: $n_attempts")

    ## The attempt count is returned as well as printed, so the comparison at the
    ## end of the session can put it beside the ABC-SMC simulation count
    return samples, n_attempts
end

Main.Notebook.abc_rejection

The components are:

# Prior
## The same priors the particle MCMC session used for these three parameters.
## The six-parameter chain also estimates D_lat, α and D_imm; the three-parameter
## chain fixes them at the values above, so it is the like-for-like comparison.
## Defined once here and read by all three samplers below: rejection draws from
## them, ABC-MCMC uses the density in its accept step, and ABC-SMC uses it in the
## importance weights
const PRIORS = [truncated(Normal(3.0, 2.0), lower = 1.0),  ## R_0
                truncated(Normal(3.0, 2.0), lower = 0.5),  ## D_inf
                Beta(2, 2)]                                ## ρ

sample_from_prior() = rand.(PRIORS)

## -Inf outside the support, so this doubles as a bounds check
log_prior(θ) = sum(logpdf(PRIORS[d], θ[d]) for d in eachindex(θ))

prior_pdf(θ) = exp(log_prior(θ))

in_prior(θ) = isfinite(log_prior(θ))

# Simulator
simulate_abc(params) = simulate_seit4l_obs(params...)

# Distance functions
distance_L1(sim, obs) = sum(abs.(sim .- obs))

## A single summary statistic rather than the whole series: the total number of
## cases over the outbreak. Cheap, and much easier to land inside a window.
distance_total(sim, obs) = abs(sum(sim) - sum(obs))

function distance_L2(sim, obs)
    return sqrt(sum((sim .- obs).^2))
end

function distance_weighted(sim, obs)
    # An observation-scaled absolute distance: closer to relative error
    weights = 1.0 ./ (obs .+ 1.0)
    return sum(weights .* abs.(sim .- obs))
end

# Observed data
observed_data = flu_tdc.obs

59-element Vector{Int64}:
  0
  1
  0
 10
  6
 32
 47
 37
 29
 11
  ⋮
  2
  0
  0
  0
  0
  0
  0
  0
  1

NoteUnderstanding the distance functions

The three distance functions have different properties:

• L1 (Manhattan): sum(|sim - obs|) - Sums absolute deviations. Robust to outliers because large errors contribute linearly. A mismatch of 20 cases counts twice as much as a mismatch of 10.

• L2 (Euclidean): sqrt(sum((sim - obs)²)) - Root of squared deviations. The squared contribution of a mismatch of 20 is four times that of 10 (400 versus 100). If it is the only mismatch, the final Euclidean distance is twice as large because of the square root.

• Weighted: sum(w × |sim - obs|) where w = 1/(obs + 1) - Weights inversely by observation size, so a mismatch of 5 when obs=2 is penalised more than when obs=50. This is an observation-scaled absolute error; it is not a Poisson discrepancy. On these data the scaling is extreme: the 16 days with no cases account for 61% of the total weight, and the days with ten cases or more for 2%.

When to use each:

• L1: Good default. Robust to occasional large mismatches.
• L2: When large deviations are particularly bad (e.g., outbreak peaks).
• Weighted: When fit quality should be relative to the observed scale, and the data have few days with no cases. The comparison of distances later in the session shows what happens when they have many.

This is the same comparison as Poisson against negative binomial, which you met yesterday. A negative log-likelihood is a distance, and each observation model fixes how steeply that distance grows with one day’s discrepancy. Take an expected count of 5 and double a discrepancy from 10 cases to 20: L1 doubles, the squared error inside L2 quadruples, the Poisson negative log-likelihood multiplies by 3, and a negative binomial with dispersion 2 multiplies by 2.2. Poisson sits between the two distances, and the negative binomial lands close to L1. Near a good fit the Poisson term is about $d^2 / 2\lambda$, which is the weighted distance above with squares in place of absolute values.

So choosing a distance and choosing an observation model are one decision, made twice: how steeply should a surprising observation count against you.

Run ABC with L1 distance:

WarningChoosing the threshold

ABC rejection sampling is sensitive to the threshold choice:

• Too strict (e.g., 100): none of the 200,000 simulations in the reference table below comes this close to the data, so it is likely to exhaust the attempt budget
• Too loose (e.g., 1000): Fast, but approximate posterior is close to the prior
• In between (e.g., 200-400): Runs quickly, but at the loose end of this band the posterior still carries a good deal of the prior

We start with a relaxed threshold (230) here. The run below accepts around 7% of what it proposes. By the guidance later in this session that is loose, and the posterior still holds a good deal of the prior. The exercise that follows asks you to tighten it.

If you get 0 samples: Increase the threshold (try 280, then 330). The sampler may be struggling to find parameter combinations that produce simulations close to the data. This can happen due to random variation - ABC is inherently stochastic.

# Start with a relaxed threshold to get some samples
threshold = 230.0  # about 7% acceptance from the prior

println("Running ABC rejection sampler...")
abc_samples, abc_attempts = abc_rejection(
    sample_from_prior,
    simulate_abc,
    distance_L1,
    observed_data,
    threshold,
    200  # Target number of samples
)

println("\nObtained $(length(abc_samples)) samples")

Running ABC rejection sampler...
Acceptance rate: 6.27%
Attempts needed: 3190

Obtained 200 samples

The sampler above answers the question for one threshold, and throws each simulation away once it has tested it. Whether a draw is accepted depends only on its distance, so if you keep the simulations instead, one batch of prior draws answers the question for every threshold at once. A stored batch like this is called a reference table, and it is how rejection ABC is usually run in practice. Building this one takes about half a minute.

## A reference table: prior draws, each simulated once, with the simulation kept
n_table = 200_000
table_θ = [sample_from_prior() for _ in 1:n_table]
table_sims = [simulate_abc(θ) for θ in table_θ]

## The L1 distance of every stored simulation from the data
table_L1 = [distance_L1(sim, observed_data) for sim in table_sims]

## Rejection ABC at threshold ε, read off the table
accepted(distances, ε) = table_θ[distances .< ε]

accepted (generic function with 1 method)

The comparisons from here on draw each posterior as a smoothed density beside the prior, and summarise it by its mean, its 90% interval and its shrinkage. Shrinkage is the share of the prior’s variance that the data have removed: 0 means the posterior is as wide as the prior, and 1 means the data have pinned the parameter down. ABC-SMC particles come with weights, so every summary is weighted, and a sample with equal weights is the special case.

"A posterior to compare: a label, an n × 3 matrix of draws, their weights and a colour."
curve(label, θ::AbstractMatrix, colour; w = ones(size(θ, 1))) =
    (; label, θ, w = w ./ sum(w), colour)
curve(label, samples::AbstractVector, colour; kwargs...) =
    curve(label, permutedims(reduce(hcat, samples)), colour; kwargs...)
curve(label, chain::DataFrame, colour; kwargs...) =
    curve(label, Matrix(chain[:, [:R_0, :D_inf, :ρ]]), colour; kwargs...)

const PAR_NAMES = ["R_0", "D_inf", "ρ"]

"Print the weighted mean, 90% interval and shrinkage of each parameter, per curve."
function summarise(curves)
    for c in curves
        columns = map(1:3) do d
            x, w = c.θ[:, d], Weights(c.w)
            lo, hi = quantile(x, w, [0.05, 0.95])
            shrinkage = 1 - var(x, w; corrected = false) / var(PRIORS[d])
            string(PAR_NAMES[d], " ", round(mean(x, w), digits = 2), " (",
                   round(lo, digits = 2), " to ", round(hi, digits = 2),
                   "), shrinkage ", round(shrinkage, digits = 2))
        end
        println(c.label, "\n    ", join(columns, ";  "))
    end
end

Main.Notebook.summarise

## One colour per method, the same in every figure, with particle MCMC in dark ink
## and a one-hue ramp for a sequence of thresholds
const SERIES = ["#2a78d6", "#eb6834", "#1baf7a"]
const BLUE_RAMP = ["#86b6ef", "#3987e5", "#1c5cab", "#0d366b"]
const INK = "#2b2b2b"
const GREY = "#9a9a9a"

const PAR_LABELS = ["R₀", "D_inf", "ρ"]
const SUPPORT = [(1.0, Inf), (0.5, Inf), (0.0, 1.0)] ## where each prior is positive
const XRANGE = [(1.0, 12.0), (0.5, 7.0), (0.0, 1.0)]

commas(n) = replace(string(n), r"(?<=\d)(?=(\d{3})+$)" => ",")

## Simulation counts in thousands once they reach five digits, to keep a key narrow
short_count(n) = n < 10_000 ? commas(n) : "$(round(Int, n / 1000))k"

## Weighted Gaussian kernel density, reflected at the edges of the prior's support
## so that a curve does not dip at a boundary. The bandwidth uses the effective
## sample size, so a few heavily weighted particles are smoothed as few
function weighted_density(x, w, xs, lo, hi)
    μ = sum(w .* x)
    h = 1.06 * sqrt(sum(w .* (x .- μ) .^ 2)) * (1 / sum(abs2, w))^(-1 / 5)
    density = zeros(length(xs))
    for (xi, wi) in zip(x, w), (k, g) in enumerate(xs)
        z = pdf(Normal(xi, h), g) + pdf(Normal(2lo - xi, h), g)
        isfinite(hi) && (z += pdf(Normal(2hi - xi, h), g))
        density[k] += wi * z
    end
    return density
end

"""
The marginal posteriors against the prior, one panel per parameter. The key sits
in the top right of the D_inf panel, where every posterior here has already
fallen away.
"""
function posterior_panels(curves)
    panels = map(1:3) do d
        lo, hi = XRANGE[d]
        xs = range(lo, hi, length = 300)
        p = plot(xlabel = PAR_LABELS[d], ylabel = d == 1 ? "Density" : "",
                 xlims = (lo, hi), legend = d == 2 ? :topright : false,
                 legendfontsize = 7, foreground_color_legend = nothing,
                 gridalpha = 0.15)
        for c in curves
            plot!(p, xs, weighted_density(c.θ[:, d], c.w, xs, SUPPORT[d]...),
                  color = c.colour, linewidth = 2, label = c.label)
        end
        plot!(p, xs, pdf.(PRIORS[d], xs), color = :black, linestyle = :dash,
              label = "prior")
    end
    return plot(panels..., layout = (1, 3), size = (900, 320),
                left_margin = 3Plots.mm, bottom_margin = 4Plots.mm)
end

Main.Notebook.posterior_panels

TipExercise: What is the threshold doing?

The next cell reads the posterior off the table at three thresholds, 400, 230 and 165, and draws each beside the prior and the particle MCMC chain that estimates the same three parameters. Predict what you will see before you run it.

1. How will the posterior at 400 compare with the prior?
2. As the threshold tightens, will the three parameters narrow at the same rate?
3. What does moving from 230 to 165 cost, and how would you decide where to stop if you had no particle MCMC answer to compare against?

thresholds_shown = [400.0, 230.0, 165.0]
for ε in thresholds_shown
    println("ε = ", ε, ": ", commas(length(accepted(table_L1, ε))), " of ",
            commas(n_table), " draws accepted (",
            round(100 * mean(table_L1 .< ε), digits = 2), "%)")
end

threshold_curves = [
    [curve("ε = $(Int(ε))", accepted(table_L1, ε), colour)
     for (ε, colour) in zip(thresholds_shown, BLUE_RAMP[2:end])];
    curve("particle MCMC", pmcmc_three, INK)
]
summarise(threshold_curves)
posterior_panels(threshold_curves)

ε = 400.0: 195,932 of 200,000 draws accepted (97.97%)
ε = 230.0: 12,868 of 200,000 draws accepted (6.43%)
ε = 165.0: 387 of 200,000 draws accepted (0.19%)
ε = 400
    R_0 3.56 (1.31 to 6.45), shrinkage -0.01;  D_inf 3.39 (0.92 to 6.39), shrinkage 0.0;  ρ 0.49 (0.13 to 0.85), shrinkage 0.03
ε = 230
    R_0 5.53 (3.41 to 7.92), shrinkage 0.24;  D_inf 1.91 (0.67 to 3.7), shrinkage 0.68;  ρ 0.57 (0.33 to 0.85), shrinkage 0.48
ε = 165
    R_0 6.17 (4.28 to 8.34), shrinkage 0.38;  D_inf 1.68 (0.62 to 2.98), shrinkage 0.79;  ρ 0.67 (0.51 to 0.84), shrinkage 0.79
particle MCMC
    R_0 7.18 (5.42 to 9.16), shrinkage 0.49;  D_inf 1.38 (0.63 to 2.34), shrinkage 0.9;  ρ 0.71 (0.64 to 0.79), shrinkage 0.96

Rejection ABC posteriors read off the reference table at three thresholds, against the prior (dashed) and the particle MCMC posterior on the same three parameters.

How far can the threshold usefully go? Run the model many times at the mean of the particle MCMC posterior, which is as good a parameter set as we have, and see how far each run lands from the data:

θ_good = [mean(pmcmc_three.R_0), mean(pmcmc_three.D_inf), mean(pmcmc_three.ρ)]
good_distances = [distance_L1(simulate_abc(θ_good), observed_data) for _ in 1:2000]

println("L1 distance at the particle MCMC mean: 5%, 50% and 95% quantiles ",
        round.(Int, quantile(good_distances, [0.05, 0.5, 0.95])))
for ε in thresholds_shown
    println("  share of runs under ", ε, ": ", round(mean(good_distances .< ε), digits = 3))
end
println("closest of all ", commas(n_table), " table draws: ", minimum(table_L1))

L1 distance at the particle MCMC mean: 5%, 50% and 95% quantiles [157, 229, 312]
  share of runs under 400.0: 0.998
  share of runs under 230.0: 0.505
  share of runs under 165.0: 0.073
closest of all 200,000 table draws: 121

TipSolution

At 400 almost every draw is accepted, and the posterior is the prior: the curves lie on the dashed line. A threshold only conditions on the data once it is tight enough to turn some parameter sets away.

No. Tightening narrows $D_\mathrm{inf}$ and $\rho$ most and $R_0$ least. From 230 to 165 the shrinkage of $\rho$ rises from about 0.5 to about 0.8, and that of $R_0$ only from 0.25 to about 0.4.

The cost is steep. About 6.5% of prior draws are accepted at 230 and about 0.2% at 165, so each accepted draw at 165 takes over thirty times as many simulations. Even at the particle MCMC mean, the middle 90% of runs land between about 155 and 310 from the data, and only about one run in twelve lands under 165. The threshold is now working inside the simulator’s own noise: most runs fail it even from the best parameters we have, so each further step down costs more simulations than the last.

There is no way to pick $\epsilon$ from first principles. In practice you lower it until either the posterior stops changing or the acceptance rate becomes unaffordable, whichever comes first. ABC-SMC, in the next section, pushes the second limit further by walking the threshold down automatically.

A summary statistic in place of the series

distance_total throws away everything except the total number of cases, one number against fifty-nine. A middle course keeps a handful of features of the curve, here the total, the day and size of the peak, and the cases in the first five days:

## Total, peak day, peak size, and cases over the first five days
epidemic_summaries(data) = [sum(data), argmax(data), maximum(data), sum(data[1:5])]

## Each summary scaled by its observed value, so that no one of them dominates
function distance_summaries(sim, obs)
    s_sim = epidemic_summaries(sim)
    s_obs = epidemic_summaries(obs)
    return sum(abs.(s_sim .- s_obs) ./ (abs.(s_obs) .+ 1.0))
end

distance_summaries (generic function with 1 method)

These distances are on different scales, so one threshold for all of them would not compare like with like. Each keeps the same share of the reference table instead, the closest 0.5% of draws, which gives every summary the same simulation budget. The total can be matched exactly: we keep only the draws whose total is the observed 312, so its threshold is zero and nothing in its posterior comes from a tolerance.

"The draws in the reference table whose distances are the closest `share` of all."
closest(distances; share = 0.005) =
    table_θ[partialsortperm(distances, 1:round(Int, share * length(distances)))]

table_total = [distance_total(sim, observed_data) for sim in table_sims]
table_summaries = [distance_summaries(sim, observed_data) for sim in table_sims]

by_series = closest(table_L1)
by_total = table_θ[table_total .== 0] ## exactly 312 cases
by_summaries = closest(table_summaries)

1000-element Vector{Vector{Float64}}:
 [5.635433790947275, 0.5593023971737829, 0.7201418629576533]
 [3.533181834663071, 1.411740466925295, 0.8441174310160298]
 [4.27050332190643, 0.9075592033583697, 0.8410303611721553]
 [4.765067696241966, 1.6452803048118168, 0.855239767201792]
 [3.934100906585165, 1.2494575796583223, 0.7793343276188517]
 [6.131349169194899, 2.544385731035058, 0.6646634484891168]
 [5.451814494704164, 2.2884207958047917, 0.8153745828813261]
 [4.954064143051553, 2.3533225127129525, 0.6887737597510087]
 [7.378973941609504, 2.764097828777778, 0.786387667905338]
 [7.481032572715655, 1.584433026133371, 0.84815698464851]
 ⋮
 [6.020586050082173, 4.092580735356665, 0.7625660909406613]
 [5.403033196472429, 2.6822643126756085, 0.8658693927637496]
 [5.847478322479997, 5.014848770169062, 0.6759061550141907]
 [5.250673000696038, 3.5890748186990735, 0.8586951703518297]
 [4.8388267377950545, 2.401768615130856, 0.7011295779467696]
 [7.549299662439237, 2.5432803810445193, 0.7131770257815042]
 [4.501135584112484, 3.552093398091084, 0.7304470140669193]
 [4.593500357488325, 1.4032607994614967, 0.6041159798465727]
 [3.9056823947238595, 1.5176373075014131, 0.7006958867512851]

TipExercise: What does the summary give up?

The next two cells draw the three posteriors beside the particle MCMC chain on the same three parameters, and then $R_0$ against $\rho$ for the full series and the total. Predict the answers before you run them.

1. The total is matched exactly, with no tolerance at all. Will its posterior agree with the full-series one?
2. Which parameter will the total tell you least about?
3. Where will the four summaries land?

summary_curves = [
    curve("full series", by_series, SERIES[1]),
    curve("total = 312 (n = $(length(by_total)))", by_total, SERIES[2]),
    curve("four summaries", by_summaries, SERIES[3]),
    curve("particle MCMC", pmcmc_three, INK),
]
summarise(summary_curves)
## Every curve here is an unweighted sample, so the plain correlation applies
for c in summary_curves
    println(rpad(c.label, 44), " cor(R_0, ρ) ", round(cor(c.θ[:, 1], c.θ[:, 3]), digits = 2))
end
posterior_panels(summary_curves)

full series
    R_0 6.05 (4.08 to 8.3), shrinkage 0.35;  D_inf 1.72 (0.66 to 3.19), shrinkage 0.79;  ρ 0.65 (0.47 to 0.83), shrinkage 0.75
total = 312 (n = 225)
    R_0 4.75 (2.96 to 7.0), shrinkage 0.37;  D_inf 3.47 (1.06 to 6.52), shrinkage 0.04;  ρ 0.79 (0.68 to 0.93), shrinkage 0.88
four summaries
    R_0 5.55 (3.63 to 7.89), shrinkage 0.33;  D_inf 2.29 (0.78 to 4.17), shrinkage 0.6;  ρ 0.75 (0.59 to 0.91), shrinkage 0.82
particle MCMC
    R_0 7.18 (5.42 to 9.16), shrinkage 0.49;  D_inf 1.38 (0.63 to 2.34), shrinkage 0.9;  ρ 0.71 (0.64 to 0.79), shrinkage 0.96
full series                                  cor(R_0, ρ) -0.1
total = 312 (n = 225)                        cor(R_0, ρ) -0.61
four summaries                               cor(R_0, ρ) -0.25
particle MCMC                                cor(R_0, ρ) -0.14

Rejection ABC posteriors under three summaries of the data. The full series and the four summaries each keep the closest 0.5% of the reference table; the total is matched exactly.

p_pair = plot(xlabel = "R₀", ylabel = "ρ", xlims = XRANGE[1], ylims = (0, 1),
              legend = :bottomright, foreground_color_legend = nothing,
              gridalpha = 0.15, size = (600, 420))
for (samples, label, colour) in ((by_total, "total only", SERIES[2]),
                                 (by_series, "full series", SERIES[1]))
    scatter!(p_pair, first.(samples), last.(samples), color = colour,
             markersize = 2.5, markerstrokewidth = 0, alpha = 0.45, label = label)
end
scatter!(p_pair, pmcmc_three.R_0[1:8:end], pmcmc_three.ρ[1:8:end], color = INK,
         markersize = 1.5, markerstrokewidth = 0, alpha = 0.35,
         label = "particle MCMC")

Draws of $R_0$ and $\rho$ accepted on the total alone and on the full series, with particle MCMC draws in dark grey.

TipWhat to look for

Matching the total exactly does not rescue it. Its posterior puts $R_0$ near 4.8 against 6.0 on the full series, and $\rho$ near 0.79 against 0.65. It barely moves $D_\mathrm{inf}$ off its prior: it removes about 4% of the prior’s variance, where the full series removes nearly 80%.

The scatter shows why. Many epidemics produce 312 cases: a small $R_0$ with a high reporting rate, a larger one with a lower rate. The total cannot tell them apart, so its posterior lies along a ridge where $R_0$ and $\rho$ trade off, with a correlation near −0.6, and the full series rules out most of that ridge. $D_\mathrm{inf}$ sets how fast the epidemic moves, and the total contains no information about timing, so it tells you almost nothing about $D_\mathrm{inf}$.

The four summaries land between the two, closer to the full series on $R_0$ and $D_\mathrm{inf}$: the peak day and the first five days put some of the timing back.

A summary is sufficient when the posterior given the summary equals the posterior given the data. The cleanest check is the full-series posterior, which shares the simulator, the priors and the fixed parameters, so the summary is the only thing that differs. That posterior is itself approximate, at a threshold well above zero, and an approximation can lose information that is in the data but it cannot invent information that is not. So what it learns about $D_\mathrm{inf}$, and the ridge it rules out, are in the data and missing from the total.

The threshold and the summary are the two approximations ABC makes. More simulations get you out of the first, but nothing gets you out of the second: once the total is all you compare, no threshold, not even zero, recovers the timing you discarded. A summary that is not sufficient puts a floor under the error, and $\epsilon \to 0$ does not lift it.

Choosing a distance

The summary sets which features of the epidemic are compared, and the distance sets how the gap between them is scored. The three distances defined earlier all compare the full series, and each keeps the closest 0.5% of the reference table, as the summaries did.

table_L2 = [distance_L2(sim, observed_data) for sim in table_sims]
table_weighted = [distance_weighted(sim, observed_data) for sim in table_sims]

distance_curves = [
    curve("L1", closest(table_L1), SERIES[1]),
    curve("L2", closest(table_L2), SERIES[2]),
    curve("weighted, 1/(obs + 1)", closest(table_weighted), SERIES[3]),
    curve("particle MCMC", pmcmc_three, INK),
]
summarise(distance_curves)
posterior_panels(distance_curves)

L1
    R_0 6.05 (4.08 to 8.3), shrinkage 0.35;  D_inf 1.72 (0.66 to 3.19), shrinkage 0.79;  ρ 0.65 (0.47 to 0.83), shrinkage 0.75
L2
    R_0 5.97 (3.95 to 8.3), shrinkage 0.32;  D_inf 1.76 (0.65 to 3.18), shrinkage 0.78;  ρ 0.71 (0.55 to 0.88), shrinkage 0.79
weighted, 1/(obs + 1)
    R_0 3.77 (1.33 to 7.11), shrinkage -0.34;  D_inf 2.32 (0.66 to 4.85), shrinkage 0.34;  ρ 0.29 (0.1 to 0.61), shrinkage 0.51
particle MCMC
    R_0 7.18 (5.42 to 9.16), shrinkage 0.49;  D_inf 1.38 (0.63 to 2.34), shrinkage 0.9;  ρ 0.71 (0.64 to 0.79), shrinkage 0.96

Rejection ABC posteriors under the three distances on the full series, each keeping the closest 0.5% of the reference table.

TipExercise: Which distance is off, and why?

One of the three posteriors pulls away from the other two and from particle MCMC. Using the weights in distance_weighted and the plot of the data at the top of the session, explain what it is fitting.

TipSolution

The weighted distance puts $\rho$ near 0.3, where the other two and particle MCMC put it near 0.7, and it leaves $R_0$ wider than its prior.

Its weights are $1/(y_t + 1)$, so a day with no cases counts 48 times as much as the peak day, with 47. Sixteen of the 59 days have no cases, and between them they account for 61% of the total weight, while the days with ten cases or more account for 2%. So the distance mostly asks whether a simulation reports nothing on the days when nothing was reported, and a low reporting rate is the easiest way to do that. L1 and L2 give almost the same posterior.

Choosing a distance is choosing an observation model, as the note on the distance functions explains. A distance that scales each day by its own count puts most of its weight on the empty days, and in this outbreak those days hold the least information.

NoteABC-MCMC, and why this session does not use it

You have built Metropolis-Hastings twice in this course: once by hand in the MCMC session, and once with a particle filter supplying the likelihood in the particle MCMC session. ABC gives you a third version of the same loop. Propose a step, simulate once from the proposed parameters, and accept if that simulation lands within ε of the data. That is ABC-MCMC, introduced by Marjoram et al. (2003). The accept step is the Metropolis-Hastings one with the likelihood ratio deleted and the distance test put in its place; the prior ratio stays, and the priors here are not uniform, so nothing cancels.

"""
ABC-MCMC. Propose, test the prior ratio, then simulate once and accept if the
simulation lands within `threshold` of the data. Returns the chain and the
acceptance rate.
"""
function abc_mcmc(θ_start, simulator, distance_fn, observed_data, threshold,
                  n_iterations; proposal_sd = [0.4, 0.3, 0.05])
    chain = Vector{Vector{Float64}}(undef, n_iterations)
    θ = copy(θ_start)
    lp = log_prior(θ)
    n_accepted = 0

    for i in 1:n_iterations
        θ_proposed = θ .+ randn(length(θ)) .* proposal_sd
        lp_proposed = log_prior(θ_proposed)

        ## The prior ratio is checked first because it is free, and a proposal
        ## that fails it needs no simulation. The proposal is symmetric, so it
        ## contributes nothing to the ratio
        if isfinite(lp_proposed) && log(rand()) < lp_proposed - lp &&
           distance_fn(simulator(θ_proposed), observed_data) < threshold
            θ, lp = θ_proposed, lp_proposed
            n_accepted += 1
        end

        chain[i] = copy(θ) ## the state after the step, whether or not it moved
    end

    return chain, n_accepted / n_iterations
end

Main.Notebook.abc_mcmc

It stalls at a tight threshold, and the stochastic simulator is the reason. Run the model at one fixed parameter set and the L1 distance still moves between runs, because each run is a different epidemic: at the particle MCMC mean, the middle 90% of runs land between about 155 and 310. At a threshold near the bottom of that range a good parameter set fails the test most of the time and passes occasionally, so the chain accepts a step when it happens to draw a lucky epidemic and then cannot leave, because leaving requires another lucky epidemic from the proposed parameters. It is stuck on a realisation, and no proposal distribution addresses that, because the proposal is not what is being rejected. The exercise in Going further puts numbers on it.

ABC-SMC: Sequential Monte Carlo ABC

Both samplers so far run out of road well before the threshold is tight. Rejection sampling wastes almost everything it draws, because it searches from the prior every time and never learns where the good parameters were. ABC-MCMC does learn that, and seizes up as a result: at a tight threshold it accepts only lucky realisations and then cannot leave them. ABC-SMC (Sequential Monte Carlo ABC) gets round both by tightening the acceptance threshold across a sequence of generations, reusing the particles that survived the previous one.

The key idea, introduced by Toni et al. (2009), is:

1. Start permissive: Begin with a loose threshold that accepts many samples
2. Gradually tighten: Reduce the threshold across generations
3. Reuse good samples: after the first generation, draw a parent from the survivors in proportion to its weight and perturb it, leaving the prior behind

The aim is to reach a threshold out of reach of rejection sampling, by refining particles that are already close. How much that saves depends entirely on the problem, the schedule and the kernel, so treat quoted speedup factors with suspicion until you have counted your own.

NoteImplementation: ABC-SMC

This is a simplified ABC-SMC for pedagogical purposes. It does compute the importance weights, which is what makes the sample valid, but production implementations (like ApproxBayes.jl or GpABC.jl) add adaptive thresholds and better perturbation kernels.

A perturbation landing outside the prior sends the loop back to its start, which draws a new parent as well as a new perturbation, so the proposal is the whole mixture restricted to the prior’s support. The constant that restriction introduces is the same for every particle, so it cancels when the weights are normalised, and the weights computed here are exact. Had the redraw kept the same parent and only re-perturbed, it would not cancel, and children of near-boundary particles would come out over-weighted.

Note on Threads.@threads: The code below uses Threads.@threads to parallelise the particle loop. The macro distributes iterations across available CPU threads, so each particle can be simulated concurrently. For this to have any effect, Julia must be started with multiple threads (e.g., julia -t auto). See the threading note in MCMC diagnostics for setup instructions. Without multiple threads, the loop runs sequentially (slower but still correct).

"Standard deviation of the perturbation kernel applied to a particle."
kernel_sd(θ) = 0.1 .* abs.(θ)

"Density of the perturbation kernel that moves `θ_from` to `θ_to`."
function kernel_pdf(θ_to, θ_from)
    sds = kernel_sd(θ_from)
    return prod(pdf(Normal(θ_from[d], sds[d]), θ_to[d]) for d in eachindex(θ_to))
end

"""
Simple ABC-SMC implementation.
Uses a sequence of decreasing thresholds, perturbing the particles that survived
the previous one. Returns one entry per generation, holding its threshold, its
particles and their importance weights, the number of simulations it cost, and
the effective sample size of its weights.
"""
function abc_smc(prior_sampler, in_prior, simulator, distance_fn, observed_data,
                thresholds, n_particles; max_simulations_per_generation = 100_000)

    particles = Vector{Vector{Float64}}()
    weights = Float64[]
    generations = []

    for (t_idx, threshold) in enumerate(thresholds)
        println("\nGeneration $t_idx: threshold = $threshold")
        new_particles = Vector{Vector{Float64}}(undef, n_particles)
        simulations = Threads.Atomic{Int}(0)

        Threads.@threads for i in 1:n_particles
            while simulations[] < max_simulations_per_generation
                θ = if t_idx == 1
                    prior_sampler()
                else
                    # Draw a survivor in proportion to its weight, then perturb it
                    θ_old = particles[sample(1:length(particles), Weights(weights))]
                    θ_old .+ randn(length(θ_old)) .* kernel_sd(θ_old)
                end

                # A perturbation can leave the prior's support; redraw rather than
                # simulate, and do not count it as a simulation
                in_prior(θ) || continue

                Threads.atomic_add!(simulations, 1)
                if distance_fn(simulator(θ), observed_data) < threshold
                    new_particles[i] = θ
                    break
                end
            end
        end

        if !all(i -> isassigned(new_particles, i), 1:n_particles)
            error("gave up at threshold $threshold after $(simulations[]) " *
                  "simulations; loosen the schedule or raise the budget")
        end

        # Importance weights. A particle reached this generation through the
        # perturbation kernel rather than from the prior, and the weights undo
        # that: numerator is the prior density, denominator is how likely the
        # previous generation was to propose it. These priors are not uniform, so
        # the numerator has to be evaluated for each particle; with a uniform
        # prior it would be a constant that divided out. Skipping the weights
        # altogether leaves the sample biased towards wherever the kernel happens
        # to concentrate.
        new_weights = if t_idx == 1
            fill(1.0 / n_particles, n_particles)
        else
            w = [prior_pdf(new_particles[i]) /
                 sum(weights[j] * kernel_pdf(new_particles[i], particles[j])
                     for j in eachindex(particles))
                 for i in 1:n_particles]
            w ./ sum(w)
        end

        particles = collect(new_particles)
        weights = new_weights
        ess = 1 / sum(abs2, weights)
        push!(generations, (; threshold, particles, weights,
                            simulations = simulations[], ess))

        weighted_R_0 = sum(weights[i] * particles[i][1] for i in eachindex(particles))
        println("  Simulations: $(simulations[]),  weighted mean R_0: ",
                round(weighted_R_0, digits=2),
                ",  ESS: ", round(ess, digits=1), " of $n_particles")
    end

    return generations
end

Main.Notebook.abc_smc

We run it twice: with 100 particles down to 165, and with 1000 particles carried on to 150, which rejection sampling reaches only at great cost.

println("Running ABC-SMC with 100 particles...")
smc_100 = abc_smc(sample_from_prior, in_prior, simulate_abc, distance_L1,
                  observed_data, [230.0, 190.0, 165.0], 100)

println("\nRunning ABC-SMC with 1000 particles...")
smc_1000 = abc_smc(sample_from_prior, in_prior, simulate_abc, distance_L1,
                   observed_data, [230.0, 200.0, 180.0, 165.0, 150.0], 1000;
                   max_simulations_per_generation = 500_000)

"The generation of an ABC-SMC run at threshold `ε`."
at_threshold(generations, ε) = generations[findfirst(g -> g.threshold == ε, generations)]

"The simulations an ABC-SMC run spent reaching `ε`, every earlier generation included."
simulations_to(generations, ε) =
    sum(g.simulations for g in generations if g.threshold >= ε)

Running ABC-SMC with 100 particles...

Generation 1: threshold = 230.0
  Simulations: 1557,  weighted mean R_0: 5.33,  ESS: 100.0 of 100

Generation 2: threshold = 190.0
  Simulations: 1170,  weighted mean R_0: 5.81,  ESS: 77.1 of 100

Generation 3: threshold = 165.0
  Simulations: 3731,  weighted mean R_0: 6.23,  ESS: 46.0 of 100

Running ABC-SMC with 1000 particles...

Generation 1: threshold = 230.0
  Simulations: 15111,  weighted mean R_0: 5.55,  ESS: 1000.0 of 1000

Generation 2: threshold = 200.0
  Simulations: 6395,  weighted mean R_0: 5.81,  ESS: 749.3 of 1000

Generation 3: threshold = 180.0
  Simulations: 11845,  weighted mean R_0: 5.96,  ESS: 668.6 of 1000

Generation 4: threshold = 165.0
  Simulations: 27398,  weighted mean R_0: 6.16,  ESS: 697.1 of 1000

Generation 5: threshold = 150.0
  Simulations: 92383,  weighted mean R_0: 6.26,  ESS: 567.2 of 1000

Main.Notebook.simulations_to

Comparing rejection ABC and ABC-SMC

Rejection sampling from the reference table and ABC-SMC at 165 aim at the same posterior: the same prior, the same distance and the same threshold. They should agree, then, and whatever separates them comes from the sampler.

rejection_165 = accepted(table_L1, 165.0)
smc_100_165 = at_threshold(smc_100, 165.0)
smc_1000_165 = at_threshold(smc_1000, 165.0)

smc_label(n, gens, g) = "SMC, $n: $(short_count(simulations_to(gens, 165.0))) sims, " *
                        "ESS $(round(Int, g.ess))"
sampler_curves = [
    curve("rejection: $(short_count(n_table)) sims, $(length(rejection_165)) draws",
          rejection_165, SERIES[1]),
    curve(smc_label(100, smc_100, smc_100_165), smc_100_165.particles, SERIES[2];
          w = smc_100_165.weights),
    curve(smc_label(1000, smc_1000, smc_1000_165), smc_1000_165.particles, SERIES[3];
          w = smc_1000_165.weights),
    curve("particle MCMC", pmcmc_three, INK),
]
summarise(sampler_curves)
posterior_panels(sampler_curves)

rejection: 200k sims, 387 draws
    R_0 6.17 (4.28 to 8.34), shrinkage 0.38;  D_inf 1.68 (0.62 to 2.98), shrinkage 0.79;  ρ 0.67 (0.51 to 0.84), shrinkage 0.79
SMC, 100: 6,458 sims, ESS 46
    R_0 6.23 (4.71 to 8.43), shrinkage 0.24;  D_inf 1.67 (0.64 to 2.86), shrinkage 0.78;  ρ 0.66 (0.5 to 0.81), shrinkage 0.8
SMC, 1000: 61k sims, ESS 697
    R_0 6.16 (4.18 to 8.48), shrinkage 0.32;  D_inf 1.75 (0.68 to 3.19), shrinkage 0.78;  ρ 0.67 (0.5 to 0.86), shrinkage 0.78
particle MCMC
    R_0 7.18 (5.42 to 9.16), shrinkage 0.49;  D_inf 1.38 (0.63 to 2.34), shrinkage 0.9;  ρ 0.71 (0.64 to 0.79), shrinkage 0.96

Rejection sampling from the reference table and two ABC-SMC runs, with 100 and 1000 particles, at the same threshold of 165. Each key entry gives the simulations spent and the effective sample size.

Rejection and the 1000-particle run draw the same curves, which is the check that the importance weights are doing their job. The 100-particle run is the rough one. The effective sample size in its label is why: the first generation is drawn from the prior and weighted equally, and each later one loses some of its effective size as the weights become uneven. That is weight degeneracy, and at 165 it leaves the 100-particle run with too few effective draws for a smooth density. Our sampler draws parents in proportion to weight, which spreads the next generation over the particles that hold it, but it never resets the weights when the effective sample size falls, and it never adapts the kernel width. Production implementations do both.

The 1000-particle run continues below 165:

ladder = [230.0, 180.0, 165.0, 150.0]
ladder_curves = [
    [curve("ε = $(Int(ε)): $(short_count(simulations_to(smc_1000, ε))) sims, " *
           "ESS $(round(Int, at_threshold(smc_1000, ε).ess))",
           at_threshold(smc_1000, ε).particles, colour;
           w = at_threshold(smc_1000, ε).weights)
     for (ε, colour) in zip(ladder, BLUE_RAMP)];
    curve("particle MCMC", pmcmc_three, INK)
]
summarise(ladder_curves)
posterior_panels(ladder_curves)

ε = 230: 15k sims, ESS 1000
    R_0 5.55 (3.37 to 7.85), shrinkage 0.26;  D_inf 1.9 (0.64 to 3.71), shrinkage 0.7;  ρ 0.57 (0.32 to 0.86), shrinkage 0.47
ε = 180: 33k sims, ESS 669
    R_0 5.96 (4.0 to 8.26), shrinkage 0.33;  D_inf 1.73 (0.64 to 3.16), shrinkage 0.78;  ρ 0.66 (0.46 to 0.87), shrinkage 0.69
ε = 165: 61k sims, ESS 697
    R_0 6.16 (4.18 to 8.48), shrinkage 0.32;  D_inf 1.75 (0.68 to 3.19), shrinkage 0.78;  ρ 0.67 (0.5 to 0.86), shrinkage 0.78
ε = 150: 153k sims, ESS 567
    R_0 6.26 (4.33 to 8.4), shrinkage 0.41;  D_inf 1.61 (0.62 to 3.01), shrinkage 0.82;  ρ 0.68 (0.54 to 0.83), shrinkage 0.84
particle MCMC
    R_0 7.18 (5.42 to 9.16), shrinkage 0.49;  D_inf 1.38 (0.63 to 2.34), shrinkage 0.9;  ρ 0.71 (0.64 to 0.79), shrinkage 0.96

The 1000-particle ABC-SMC run at four of its thresholds. Each key entry gives the simulations spent reaching that threshold and the effective sample size there.

How much ABC-SMC saves depends on how rare acceptance is at the threshold you want. The reference table measures that directly, as the share of prior draws that land under each threshold.

grid = 140:5:400
acceptance = [mean(table_L1 .< ε) for ε in grid]

plot(grid, 100 .* max.(acceptance, 1e-4), yscale = :log10,
     xlabel = "Threshold ε", ylabel = "Prior acceptance (%)",
     label = "", linewidth = 2)
vline!([230.0], label = "first rejection run", linestyle = :dash)
vline!([165.0, 150.0], label = "ABC-SMC comparisons", linestyle = :dot)

The share of prior draws in the reference table that land under each threshold, on a log scale.

Acceptance falls away far faster than the threshold does, so the cell below turns the curve into the number that matters: how many simulations each method spends per effective sample.

for ε in (230.0, 165.0, 150.0)
    prior_acceptance = mean(table_L1 .< ε)

    ## Rejection draws are unweighted, so its effective sample size is its draw count
    rejection_per_ess = 1 / prior_acceptance

    ## The ABC-SMC particles are weighted, and the earlier generations cannot be
    ## skipped, so all of them count against the effective sample size at ε
    smc_per_ess = simulations_to(smc_1000, ε) / at_threshold(smc_1000, ε).ess

    println("At ε = ", ε, ", simulations per effective sample: rejection ",
            commas(round(Int, rejection_per_ess)), " (prior acceptance ",
            round(100 * prior_acceptance, digits = 3), "%), ABC-SMC ",
            commas(round(Int, smc_per_ess)))
end

At ε = 230.0, simulations per effective sample: rejection 16 (prior acceptance 6.434%), ABC-SMC 15
At ε = 165.0, simulations per effective sample: rejection 517 (prior acceptance 0.194%), ABC-SMC 87
At ε = 150.0, simulations per effective sample: rejection 1,905 (prior acceptance 0.052%), ABC-SMC 270

At 230 the two cost the same, about 15 simulations per effective draw, because 6.5% of prior draws are accepted there and ABC-SMC’s first generation is rejection sampling from the prior. At 165 ABC-SMC needs about 90 simulations per effective draw against about 500 for rejection, and at 150 about 250 against about 2,100, so its advantage grows as the threshold tightens. The factor moves between runs, so take the ordering and the reason for it: ABC-SMC concentrates its effort where the data are without being told where that is, and that wins whenever prior acceptance at the threshold you want is small.

Two things flatter ABC-SMC in this count. Its effective sample size comes from the weights alone, and ignores that particles descended from the same parent sit close together. And the reference table can be read at any threshold, summary or distance, where ABC-SMC has to run again for each, so at a threshold the prior still reaches, rejection sampling from a table is the better tool.

ImportantAssessing ABC posterior quality

Rejection ABC does not produce chains, so chain diagnostics like R-hat have nothing to work with. ABC-SMC produces weights, and the effective sample size of those weights tracks the same problem, as above. Beyond that, assess ABC posterior quality by:

• Visual comparison with likelihood-based posteriors: When a likelihood is available (as in this session), run both methods on the same parameters and compare. Large discrepancies suggest either the ABC threshold is too loose or the distance metric captures different features than the likelihood. With different parameters fixed, the two are aiming at different posteriors, and the comparison below shows how far apart that alone can put them.
• Threshold sensitivity: Run ABC with several thresholds. If the posterior changes substantially as you tighten the threshold, it has not yet converged, and you need a stricter threshold and more computation.
• Shrinkage against the prior: A parameter whose shrinkage is near zero is not being informed by the data, whether because the threshold is too loose or because the summary discards what would inform it.
• Replicate runs: Run ABC multiple times with different random seeds. If posteriors differ substantially between runs, you need more particles or a less strict threshold.

Do they give similar results?

The ABC posterior we take forward is the last ABC-SMC generation, 1000 particles at a threshold of 150, the tightest in this session. The figure holds it against both particle MCMC chains: the one on the same three parameters, and the one from the previous session on all six.

abc_last = smc_1000[end]
final_curves = [
    curve("ABC-SMC, ε = $(Int(abc_last.threshold))", abc_last.particles, SERIES[1];
          w = abc_last.weights),
    curve("particle MCMC, 3 parameters", pmcmc_three, INK),
    curve("particle MCMC, 6 parameters", pmcmc_chain, GREY),
]
summarise(final_curves)
posterior_panels(final_curves)

ABC-SMC, ε = 150
    R_0 6.26 (4.33 to 8.4), shrinkage 0.41;  D_inf 1.61 (0.62 to 3.01), shrinkage 0.82;  ρ 0.68 (0.54 to 0.83), shrinkage 0.84
particle MCMC, 3 parameters
    R_0 7.18 (5.42 to 9.16), shrinkage 0.49;  D_inf 1.38 (0.63 to 2.34), shrinkage 0.9;  ρ 0.71 (0.64 to 0.79), shrinkage 0.96
particle MCMC, 6 parameters
    R_0 6.32 (4.53 to 8.34), shrinkage 0.46;  D_inf 2.12 (0.82 to 3.55), shrinkage 0.76;  ρ 0.69 (0.61 to 0.77), shrinkage 0.95

The final ABC-SMC generation against particle MCMC on the same three parameters and on all six.

The exercise below asks about the latent period, which the six-parameter chain estimates and the other two fix, so here it is beside the sum of the latent and infectious periods.

D_lat_six = pmcmc_chain.D_lat
println("The six-parameter chain puts D_lat at ", round(mean(D_lat_six), digits = 2),
        " (", round(quantile(D_lat_six, 0.05), digits = 2), " to ",
        round(quantile(D_lat_six, 0.95), digits = 2), "), where the other two fix it at ",
        D_LAT)

println("D_lat + D_inf: six-parameter chain ",
        round(mean(D_lat_six .+ pmcmc_chain.D_inf), digits = 2),
        ", three-parameter chain ", round(D_LAT + mean(pmcmc_three.D_inf), digits = 2))

The six-parameter chain puts D_lat at 1.34 (0.66 to 2.04), where the other two fix it at 2.0
D_lat + D_inf: six-parameter chain 3.45, three-parameter chain 3.38

The two particle MCMC chains separate two things that a single comparison runs together. Where ABC differs from the chain on the same three parameters, the difference is the approximation: the threshold and the distance. Where the two chains differ from each other, the difference is the choice of which parameters to fix, and ABC has no part in it.

TipExercise: Where do they disagree, and why?

1. Compare ABC with the particle MCMC chain on the same three parameters. Where do they disagree, in the centre of a posterior or in its width?
2. Compare the two particle MCMC chains. Which parameter moves most when $D_\mathrm{lat}$, $\alpha$ and $D_\mathrm{imm}$ are estimated, and why that one?
3. Had you compared ABC with the six-parameter chain alone, what would you have concluded about its $D_\mathrm{inf}$?

TipSolution

ABC sits below the chain on the same three parameters for $R_0$, 6.3 against 7.2, and every ABC interval is wider, most of all that of $\rho$, which is about twice as wide. That is the approximation. At a threshold of 150 ABC still keeps simulations that miss the data by a wide margin, and the closest of the 200,000 table draws came no nearer than 121. The L1 distance decides which of those it keeps, where particle MCMC weighs every path by the Poisson likelihood, and $R_0$ has nearly stopped moving between 165 and 150, so that gap will not close cheaply.

The two particle MCMC chains differ most on $D_\mathrm{inf}$: 2.12 when the latent period is estimated, and 1.38 when it is fixed at 2. The six-parameter chain puts $D_\mathrm{lat}$ at 1.34, and what the daily counts constrain is how quickly the epidemic moves, which depends on the two periods together. Their sum is about 3.4 in both chains, so fixing $D_\mathrm{lat}$ at 2, longer than the data prefer, pulls $D_\mathrm{inf}$ down to keep the sum, and $R_0$ moves up with it. $\rho$ barely moves.

Against the six-parameter chain alone, ABC’s $D_\mathrm{inf}$ of 1.6 looks too short, and ABC would take the blame for a gap that fixing $D_\mathrm{lat}$ created. Against the chain with the same target it is slightly longer, by about 0.2. An approximate posterior needs a reference with the same target before you can tell its errors apart from your modelling choices.

Learning points

TipLearning points

• ABC replaces the likelihood with an explicit distance, but choosing L1 over L2 over a weighted distance is the same kind of choice as picking Poisson over negative binomial. A distance that scales each day by its own count put most of its weight on the days with no cases, and pulled the reporting rate down to fit them
• The acceptance threshold is the hard part: too tight and nothing is accepted, too loose and you are sampling the prior. A reference table of stored simulations gives rejection ABC at every threshold, summary and distance from one batch
• ABC-SMC walks the threshold down in stages, concentrating effort where the data are without being told where that is, with weights that become uneven. It wins where prior acceptance at the threshold you want is small, and earns nothing at a loose threshold the prior already reaches
• Comparing a summary instead of the whole series is the concession you cannot undo. A tighter threshold recovers what the tolerance gave up; nothing recovers what the summary discarded
• Hold an approximate posterior against a reference with the same target, or you cannot tell its errors from your modelling choices
• The particle filter fits this model better than ABC does. ABC earns its place where the filter cannot run: no observation density to weight with, or a simulator you cannot stop and restart

Going further

ABC-MCMC

The callout earlier in this session defines abc_mcmc. Run it down a ladder of thresholds, and then sweep the proposal width at the loosest and the tightest, reporting the effective sample size of the $R_0$ draws alongside the acceptance rate.

TipExercise: why tuning will not save this

In the MCMC session a low acceptance rate meant the proposal was too wide, and you fixed it by tuning proposal_sd. Predict what the sweep will show at ε = 230 and at ε = 165 before you run the cells, then explain the difference between the two columns.

start_distances = [distance_L1(simulate_abc(θ), observed_data) for θ in abc_samples]
θ_start = abc_samples[argmin(start_distances)]

"Effective sample size of the R_0 draws in an ABC-MCMC chain."
ess_R_0(chain) = ess(reshape([θ[1] for θ in chain], :, 1)) ## draws by chains

for ε in (230.0, 200.0, 180.0, 165.0)
    chain, acceptance = abc_mcmc(θ_start, simulate_abc, distance_L1,
                                 observed_data, ε, 4000)
    println("ε = ", ε,
            ": acceptance ", lpad(round(100 * acceptance, digits = 1), 5),
            "%, ESS ", lpad(round(Int, ess_R_0(chain)), 3), " from 4000 draws")
end

ε = 230.0: acceptance  34.1%, ESS  12 from 4000 draws
ε = 200.0: acceptance  21.6%, ESS  22 from 4000 draws
ε = 180.0: acceptance   9.0%, ESS  10 from 4000 draws
ε = 165.0: acceptance   3.2%, ESS  14 from 4000 draws

for scale in (0.25, 1.0, 2.5, 5.0)
    proposal_sd = scale .* [0.4, 0.3, 0.05]
    for ε in (230.0, 165.0)
        chain, acceptance = abc_mcmc(θ_start, simulate_abc, distance_L1,
                                     observed_data, ε, 4000;
                                     proposal_sd = proposal_sd)
        println("proposal × ", lpad(scale, 4), ", ε = ", ε,
                ": acceptance ", lpad(round(100 * acceptance, digits = 1), 5),
                "%, ESS ", lpad(round(Int, ess_R_0(chain)), 3))
    end
end

proposal × 0.25, ε = 230.0: acceptance  41.8%, ESS  10
proposal × 0.25, ε = 165.0: acceptance   3.3%, ESS  10
proposal ×  1.0, ε = 230.0: acceptance  39.2%, ESS  25
proposal ×  1.0, ε = 165.0: acceptance   3.4%, ESS  10
proposal ×  2.5, ε = 230.0: acceptance  21.4%, ESS  41
proposal ×  2.5, ε = 165.0: acceptance   2.1%, ESS  13
proposal ×  5.0, ε = 230.0: acceptance   8.2%, ESS  88
proposal ×  5.0, ε = 165.0: acceptance   0.5%, ESS  30

At ε = 230 the effective sample size rises as the proposal widens and then falls again once too little is accepted, so there is an optimum you can find. At ε = 165 there is nothing to find: the acceptance rate falls and the effective sample size sits between about ten and twenty out of four thousand whatever you set.

Summary statistics for ABC

The four summaries in the main session kept some of the timing and lost the rest. Others worth trying on epidemic data include the initial growth rate, the day the epidemic ends and the final attack rate. Choosing summaries by hand is problem-specific and takes experimentation, so a common alternative is to learn them: regress the parameters on the simulated data in a pilot run, and use the fitted values as the summaries. That is semi-automatic ABC (Fearnhead and Prangle, 2012), and the reference table built earlier is exactly the pilot run it needs.

When the simulator is noisy

Every distance in this session had to clear the noise in the simulator: even at the parameters the particle filter favours, most runs land well away from the data. That forces a loose threshold, and a loose threshold pulls the posterior back towards the prior. Several methods take this on directly:

• Regression adjustment (Beaumont, Zhang and Balding, 2002) accepts at a loose threshold, then corrects each accepted parameter set by regressing the parameters on how far its summaries missed the data.
• Several simulations per parameter set help ABC-MCMC: use the share of runs that land within the threshold as the likelihood in the accept step. The estimate is unbiased, so the target does not change and the chain sticks less, which is the same construction as particle MCMC.
• Synthetic likelihood (Wood, 2010) simulates many times at each parameter set, fits a Gaussian to the summaries, and uses its density as the likelihood.
• Neural posterior estimation, the subject of its own session, learns the posterior as a function of the data from simulations.
• Conditioning as the epidemic runs is what the particle filter does: it drops simulations that drift away from the data day by day, where ABC judges each run whole at the end. That is much of why particle MCMC fits this outbreak better.

Robust distance metrics

For ABC with outliers or model misspecification, consider robust distances:

• Median absolute deviation in place of mean squared error
• Quantile-based distances that down-weight extreme values
• Feature-based distances that focus on epidemiologically meaningful patterns

Try replacing distance_L1 in abc_rejection with a median-based version, distance_median(sim, obs) = median(abs.(sim .- obs)), and compare the posteriors. A robust distance stops a handful of large daily discrepancies dominating the fit, so how much the posterior moves tells you how much of it was being decided by a few days.

Next session

You have now met every family of method the course covers. Synthesis and next steps puts them side by side, asks which one your own problem needs, and points at the further topics that go beyond it.

References

• Beaumont, M.A., Zhang, W. and Balding, D.J. (2002). Approximate Bayesian computation in population genetics. Genetics 162(4), 2025-2035. (Regression adjustment)
• Marjoram, P. et al. (2003). Markov chain Monte Carlo without likelihoods. PNAS 100(26), 15324-15328. (ABC-MCMC)
• Toni, T. et al. (2009). Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J. R. Soc. Interface 6(31), 187-202. (The ABC-SMC paper - highly readable introduction)
• Wood, S.N. (2010). Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466, 1102-1104. (Synthetic likelihood approach)
• Marin, J.-M. et al. (2012). Approximate Bayesian computational methods. Statistics and Computing 22, 1167-1180.
• Fearnhead, P. and Prangle, D. (2012). Constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation. Journal of the Royal Statistical Society: Series B 74(3), 419-474.
• Cranmer, K. et al. (2020). The frontier of simulation-based inference. PNAS 117(48), 30055-30062. (Comprehensive review of modern likelihood-free methods)