Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) is what you reach for when the likelihood $p(y\mid \theta )$ is intractable (or you simply don't want to write it down), but you can simulate synthetic data from the model. ABC replaces likelihood evaluation with simulate-and-compare: draw $\theta$ from the prior, simulate $y_{\text{sim}}$, and accept $\theta$ if $y_{\text{sim}}$ is close enough to the real observation $y_{\text{obs}}$.

The model: simulate, don't observe

/// `mu ~ Normal(0, 1)`; forward-simulate `y_sim ~ Normal(mu, 0.5)` (a `sample`,
/// not an `observe` -- ABC never scores this density, only compares outcomes).
fn model() -> Model<f64> {
    sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
        .bind(|mu| sample(addr!("y_sim"), Normal::new(mu, 0.5).unwrap()).map(move |_| mu))
}

Notice the second line uses sample, not observe: ABC never scores a likelihood, so the model's job is to forward-simulate a synthetic observation, not condition on a real one. The address y_sim is read back out of the trace after each prior draw.

Why this converges to the same posterior

Accepting a draw when $|y_{\text{sim}}-y_{\text{obs}}|\le \epsilon$ and shrinking $\epsilon \to 0$ converges to conditioning on $y_{\text{sim}}=y_{\text{obs}}$ exactly — the same event a direct Bayesian update on $y\sim \mathcal{N}(\mu ,0.5^2)=1.5$ conditions on. So in the small-tolerance limit, ABC targets the same $\mathcal{N}(1.2,0.2)$ posterior as the SMC and VI tutorials. At any finite tolerance, though, ABC is only an approximation of that target — the price paid for not needing a likelihood at all.

Running ABC-SMC

A single-shot rejection ABC (abc_rejection) at a small tolerance can require an enormous number of prior draws to find even one acceptance. abc_smc_weighted fixes this the same way adaptive_smc does: a schedule of shrinking tolerances, with each stage's population built by perturbing and re-weighting the previous one (Beaumont et al. 2009 / Toni et al. 2009) rather than restarting from the prior every time.

    let observed: Vec<f64> = vec![1.5];
    let mut rng = StdRng::seed_from_u64(7);

    let config = ABCSMCConfig {
        initial_tolerance: 2.0,
        // Geometrically shrinking tolerance schedule (Beaumont/Toni-style):
        // each stage's population is built by perturbing and re-weighting the
        // previous one (see `abc_smc_weighted`'s rustdoc).
        tolerance_schedule: vec![1.0, 0.5, 0.25, 0.1],
        particles_per_round: 500,
    };

    let result = abc_smc_weighted(
        &mut rng,
        model,
        // The simulator just reads back the model's own forward-simulated
        // synthetic observation.
        |trace| vec![trace.get_f64(&addr!("y_sim")).unwrap()],
        &observed,
        &EuclideanDistance,
        config,
        // Attempt budget per stage (finding FG-34: bounded, typed-error
        // instead of an unbounded loop / panic on an empty population).
        200_000,
    )
    .expect("ABC-SMC should complete with this many particles/attempts");

Each stage is bounded by an attempt budget (200_000 here): if a stage can't fill its particle quota within that budget, abc_smc_weighted returns a typed ABCError (EmptyInitialPopulation or StageExhausted) instead of looping forever or panicking on an empty population.

Reading the result

abc_smc_weighted returns an ABCSMCResult with a weighted particle population (unlike abc_smc, its equally-weighted convenience wrapper) plus a weighted_mean helper:

    let posterior_mean = result
        .weighted_mean(&addr!("mu"))
        .expect("mu is present in every particle's trace");
    let posterior_var: f64 = {
        let mut num = 0.0;
        let mut den = 0.0;
        for p in &result.particles {
            let mu = p.trace.get_f64(&addr!("mu")).unwrap();
            num += p.weight * (mu - posterior_mean).powi(2);
            den += p.weight;
        }
        num / den
    };
    println!("Posterior mean(mu) ~= {posterior_mean:.4} (target: {POSTERIOR_MEAN})");
    println!(
        "Posterior sd(mu)   ~= {:.4} (target: {POSTERIOR_SD:.4})",
        posterior_var.sqrt()
    );

Checking against ground truth

    // ABC is only asymptotically (tolerance -> 0) exact, so these bands are
    // deliberately wider than the SMC example's: at final tolerance 0.1 with
    // 500 particles, the Beaumont-kernel ABC-SMC posterior mean should still
    // land within a few tenths of the exact value, but a broken importance
    // correction (e.g. reverting to the prior-replacement heuristic of
    // finding FG-09) biases it by much more than this band allows.
    assert!(
        (posterior_mean - POSTERIOR_MEAN).abs() < 0.3,
        "ABC posterior mean {posterior_mean} too far from target {POSTERIOR_MEAN}"
    );
    assert!(
        posterior_var.sqrt() < POSTERIOR_SD * 2.0,
        "ABC posterior sd {} implausibly larger than target {POSTERIOR_SD}",
        posterior_var.sqrt()
    );

Note the wider tolerance band compared to the SMC tutorial's assertions — that gap is the ABC trade-off made visible: approximate inference, in exchange for applicability to simulators with no tractable likelihood whatsoever.

Run it yourself:

cargo run --example abc_inference

Next

• Variational Inference — same target again, this time fit by optimization.
• Back to Sequential Monte Carlo if you haven't read it yet.