Sequential Monte Carlo

Sequential Monte Carlo (SMC) maintains a population of weighted particles and moves them through a sequence of intermediate target distributions, resampling and rejuvenating along the way, until the population approximates the posterior. Unlike MCMC's single evolving chain, SMC gives you many (weakly correlated) draws per run and an unbiased estimate of the log marginal likelihood — useful for model comparison, which no single-chain MCMC method provides directly.

The model

/// `mu ~ Normal(0, 1)`; observe `y ~ Normal(mu, 0.5)` at the fixed value 1.5.
fn model() -> Model<f64> {
    sample(addr!("mu"), Normal::new(0.0, 1.0).unwrap())
        .bind(|mu| observe(addr!("y"), Normal::new(mu, 0.5).unwrap(), 1.5).map(move |_| mu))
}

This is the conjugate Normal-Normal setup introduced in the section overview: prior $\mu \sim \mathcal{N}(0,1)$, likelihood $y\mid \mu \sim \mathcal{N}(\mu ,0.5^2)$ observed at $y=1.5$, with closed-form posterior $\mathcal{N}(1.2,0.2)$.

How adaptive_smc gets there

adaptive_smc targets the sequence of likelihood-tempered distributions

$${\pi }_{\beta }(\theta )\propto p(\theta )\cdot p(y\mid \theta {)}^{\beta },\beta :0\to 1$$

so ${\pi }_0$ is the prior (trivial to sample) and ${\pi }_1$ is the posterior. At each step it:

1. Adapts the next $\beta$ by bisection so the reweighted effective sample size (ESS) hits the configured threshold — a big jump when particles agree, a small one when they don't.
2. Reweights particles by the incremental likelihood factor and folds the reweighting into a running log-evidence estimate.
3. Resamples (when not the final step) to discard low-weight particles, then rejuvenates with a handful of ${\pi }_{\beta }$-invariant Metropolis-Hastings moves to restore diversity lost during resampling.

    let mut rng = StdRng::seed_from_u64(42);
    let config = SMCConfig {
        resampling_method: ResamplingMethod::Systematic,
        ess_threshold: 0.5,
        // Rejuvenation moves diversify particles after each resample, which is
        // what makes this genuine (multi-step) tempered SMC rather than a
        // single importance-sampling reweight (see `adaptive_smc`'s rustdoc).
        rejuvenation_steps: 3,
    };
    let num_particles = 2000;
    let result = adaptive_smc(&mut rng, num_particles, model, config);

rejuvenation_steps: 3 is what makes this genuine multi-step tempered SMC rather than a single importance-sampling reweight — with zero rejuvenation steps, particle positions never change between resamples, so adaptive_smc degenerates to one prior-to-posterior importance-sampling jump (still correct, just less effective for hard posteriors).

Reading the result

adaptive_smc returns an SMCResult, which dereferences to Vec<Particle> (so slice/iterator methods work directly) and additionally carries log_evidence:

    // Weighted posterior mean/variance over `mu` from the final particle population.
    let weighted_mean: f64 = result
        .iter()
        .filter_map(|p| p.trace.get_f64(&addr!("mu")).map(|mu| p.weight * mu))
        .sum();
    let weighted_var: f64 = result
        .iter()
        .filter_map(|p| {
            p.trace
                .get_f64(&addr!("mu"))
                .map(|mu| p.weight * (mu - weighted_mean).powi(2))
        })
        .sum();
    println!("Posterior mean(mu) ~= {weighted_mean:.4} (exact: {POSTERIOR_MEAN})");
    println!("Posterior var(mu)  ~= {weighted_var:.4} (exact: {POSTERIOR_VAR})");

Checking against ground truth

    // With 2000 particles and 3 rejuvenation moves per intermediate temper
    // step, the weighted-mean Monte Carlo error is well under 0.1 for this
    // 1-D conjugate model; 0.15 gives comfortable headroom without being
    // vacuous (a broken importance weight -- e.g. reintroducing the
    // prior-squaring bug of FG-03 -- would shift the mean by several tenths).
    assert!(
        (weighted_mean - POSTERIOR_MEAN).abs() < 0.15,
        "SMC posterior mean {weighted_mean} too far from exact {POSTERIOR_MEAN}"
    );
    assert!(
        (weighted_var - POSTERIOR_VAR).abs() < 0.1,
        "SMC posterior var {weighted_var} too far from exact {POSTERIOR_VAR}"
    );
    assert!(
        result.log_evidence.is_finite(),
        "log-evidence estimate must be finite"
    );

Run it yourself:

cargo run --example smc_inference

Next

• Approximate Bayesian Computation — same target, no likelihood required at all.
• Variational Inference — same target, fit by optimization instead of sampling.