When one proposal isn't enough: march a weighted swarm of particles through a ladder of distributions, resampling to survive, and read off the normalizing constant nobody else can.
You want a single number: the probability that a standard normal random variable exceeds 4. Call it P(X > 4). You fire 10,000 independent samples at it and count how many land past 4. Easy, right?
The true value is about 3.2×10−5. So across 10,000 draws you expect 10000 × 0.000032 ≈ 0.32 hits — less than one. The most likely outcome is zero samples past 4, which makes your estimate exactly 0, with infinite relative error.
Cranking N up to a million barely helps: to even see one hit you need roughly 1/p ≈ 31,600 samples, and to get a stable estimate you need around 100/p ≈ 3 million. For X > 5 the cost is another factor of ten worse. The price of a rare event explodes as 1/p.
You met one cure in the last lesson. Importance sampling (see that lesson) tilts a single proposal out toward the tail so samples land where you care. But you must know where to aim, and in high dimensions a single tilted proposal collapses — one weight eats everything.
Here is the new idea. Instead of one giant leap to the tail, slide a whole swarm of particles outward in stages: first past 2, then past 3, then past 4. At each stage you reweight, then breed the survivors — clone the particles that made it, kill the ones that didn't. The swarm physically crawls into the tail and stays populated the whole way.
That staged march is Sequential Monte Carlo (SMC). It converts one impossible rare-event problem into a sequence of easy ones whose product is the answer. This lesson builds the static-target sampler version: tempering, resampling, evidence. (For the time-series cousin that tracks a moving hidden state, see particle filter — same engine, opposite purpose, Chapter 11.)
The faint curve is the N(0,1) density. A threshold bar sits at L. In Naive mode, dots are i.i.d. normal — count how few ever land right of the bar. In SMC mode, particles right of the bar get cloned and the rest die, then jitter, so the swarm crawls into the tail. Press Play to march L from 2 → 3 → 4.
No randomness needed — we predict the failure analytically. Let Φ be the standard normal CDF, so the tail probability is p = 1 − Φ(4) = 3.167×10−5.
The expected number of tail hits in N = 10,000 draws is N·p = 10000 × 3.167×10−5 = 0.317 — less than one hit on average, so the modal estimate is 0/10000 = 0.
Now stage the same probability as a product of three easy conditional factors. Each factor is a probability the swarm can estimate with a well-populated population:
Multiply them: 0.02275 × 0.05934 × 0.02346. The first product 0.02275 × 0.05934 = 0.0013500; times 0.02346 = 3.167×10−5 — exactly the tail probability, reconstructed. The telescoping is exact because P(X>2)·P(X>3|X>2)·P(X>4|X>3) = P(X>4) by the chain rule of conditional probability.
# Naive MC fails on the tail; staged conditional factors do not import numpy as np rng = np.random.default_rng(0) x = rng.standard_normal(10_000) print((x > 4).mean()) # almost always 0.0 from scipy.stats import norm levels = [2, 3, 4]; prev = 0; prob = 1.0 for L in levels: prob *= (1 - norm.cdf(L)) / (1 - norm.cdf(prev) if prev else 1.0) prev = L print(prob) # ~3.167e-5, reconstructs the tail
# library one-liner: the exact survival function from scipy.stats import norm; norm.sf(4) # SMC estimates each conditional factor instead
You already built the core trick in importance sampling: to estimate an expectation under a hard target π, draw from an easy proposal q, then correct each sample by its importance weight w = π(x)/q(x). That lesson owns the identity; here we only recall it, then go deep on what happens when you run many weighted particles at once.
One sample, one weight, one correction. The weight is the ratio of how likely the sample is under the target versus the proposal. If the proposal over-visits a region, weights there are small (down-weighted). If the proposal under-visits a region the target loves, weights there are large (up-weighted). The weights repair the mismatch.
But a single particle is fragile. If its weight is huge it dominates the average, and your "estimate" is essentially one number. SMC's entire premise is to run a population of weighted particles — and, crucially, to let good ones breed and bad ones die. Before we can breed them, we must be crystal-clear on what one particle's weight is and what it means.
The self-normalized estimate of an expectation Eπ[X] is the weighted average Σi wi xi / Σi wi. The denominator normalizes the weights so they sum to one — which is why any constant scaling of all weights cancels and leaves the estimate unchanged.
That self-normalization is also why a single weight is meaningless on its own. A weight of 0.22 says nothing until you compare it to the others. What matters is the normalized weight Wi = wi / Σj wj. Two runs whose weights differ by a factor of 1000 give identical estimates.
The absolute weight scale carries real information in exactly one place: when you chase the normalizing constant (Chapters 6–7). There the un-normalized weight mass is the prize. For now, for estimating expectations, only the normalized weights matter.
Target π = N(0,1) solid; proposal q = N(3,1) dashed (deliberately mismatched). Click anywhere to drop a particle at that x: vertical lines hit both curves, and the weight w = π(x)/q(x) is drawn as a circle whose radius scales with w. Particles near 0 get huge weights; particles near 3 get tiny ones. Drag the proposal-mean slider and watch all weights rescale.
Target π = N(0,1), proposal q = N(3,1). A particle lands at x = 2. The Gaussian density is φ(x; μ) = (1/√(2π)) exp(−(x−μ)2/2).
Numerator: π(2) = (1/√(2π)) exp(−22/2) = (1/√(2π)) exp(−2) = 0.05399.
Denominator: q(2) = (1/√(2π)) exp(−(2−3)2/2) = (1/√(2π)) exp(−0.5) = 0.24197.
Weight: w = π(2)/q(2) = 0.05399 / 0.24197 = 0.2231. The 1/√(2π) prefactors cancel in the ratio, so w = exp(−2)/exp(−0.5) = exp(−1.5) = 0.2231 — same number, by hand.
Interpretation: x = 2 is fairly likely under the proposal (near its mean 3) but only moderately likely under the target, so this sample is down-weighted to 0.22. The proposal over-represented this region; the weight corrects it.
# self-normalized IS estimate of E_pi[X] with a swarm of weights import numpy as np from scipy.stats import norm rng = np.random.default_rng(0) x = rng.normal(3, 1, size=1000) # draw from proposal q=N(3,1) w = norm.pdf(x, 0, 1) / norm.pdf(x, 3, 1) # weights pi/q W = w / w.sum() # normalize print((W * x).sum()) # weighted estimate of E_pi[X] ~ 0
# library one-liner W = target.pdf(x) / proposal.pdf(x); est = np.average(x, weights=W)
You start with 4 particles, all equally weighted at 0.25 — a healthy swarm. Then you observe some data and multiply each weight by how well that particle explains it. The incremental likelihood factors come out to (2.0, 0.1, 0.05, 0.01). Renormalize.
Now particle 1 holds 92.6% of the weight; the other three together hold 7.4%. Do this a few more times and one particle holds 99.99%. Your "1000-particle" estimate is secretly a one-particle estimate.
This collapse is weight degeneracy: as you sequentially multiply weights by more and more factors, the variance of the weights explodes and almost all the normalized mass collects on a single particle. It is the disease that kills naive sequential importance sampling, and the rest of this lesson is the cure.
We measure the damage with the effective sample size, ESS = 1 / Σi Wi2, where W are the normalized weights. If all N weights are equal, ESS = N (full health). If one weight is 1 and the rest are 0, ESS = 1 (total collapse). ESS answers: "how many independent samples is my weighted swarm actually worth?"
A subtle but important point: degeneracy is about the weights, not the particle locations. The particles can be beautifully spread over the target, yet ESS = 1 because one of them happened to win a slightly higher weight at every step and the ratios compounded multiplicatively.
That diagnosis is what tells us the fix can never be "move the particles." The fix is resampling (next chapter): discard the near-zero-weight particles and clone the survivors so the weights reset to uniform and ESS jumps back to N.
Bars start equal. Each Observe multiplies the weights by a fresh incremental-likelihood vector and renormalizes; the bars redistribute until one towers over the rest. The ESS meter falls from N toward 1. The slider controls how informative (spread out) each step's likelihood is — more informative steps crash ESS faster.
Take the fixed, already-skewed normalized weight vector W = (0.7, 0.1, 0.1, 0.05, 0.05). Then ESS = 1 / Σ Wi2.
Sum of squares: 0.72 + 0.12 + 0.12 + 0.052 + 0.052 = 0.49 + 0.01 + 0.01 + 0.0025 + 0.0025 = 0.515.
So ESS = 1 / 0.515 = 1.942. Five particles carry only about 1.9 particles' worth of information — the swarm is already three-fifths dead.
Now watch one observation destroy a healthy swarm. Start uniform: W = (0.25, 0.25, 0.25, 0.25), so ESS = 1/(4 × 0.0625) = 1/0.25 = 4.0. Multiply by incremental likelihood g = (2, 0.1, 0.05, 0.01): the un-normalized products are (0.5, 0.025, 0.0125, 0.0025), summing to 0.54.
Normalize: W = (0.5, 0.025, 0.0125, 0.0025) / 0.54 = (0.9259, 0.0463, 0.0231, 0.0046). New sum of squares = 0.92592 + 0.04632 + 0.02312 + 0.00462 = 0.8573 + 0.0021 + 0.0005 + 0.00002 = 0.86. ESS = 1/0.86 = 1.1627. One informative observation nearly destroyed the swarm.
# effective sample size and the collapse import numpy as np def ess(W): W = W / W.sum() return 1.0 / np.sum(W**2) # 1 / sum of squared normalized weights W = np.ones(4) / 4 # healthy: ESS=4.0 g = np.array([2, 0.1, 0.05, 0.01]) # incremental likelihood W = W * g # reweight print(np.round(W / W.sum(), 4), 'ESS=', round(ess(W), 3)) # ~1.163
# library one-liner: the trigger for resampling ess = 1.0 / np.sum((w / w.sum())**2) # resample when this drops below N/2
Degeneracy hands all the weight to one particle. The fix is brutal and biological: hold a lottery where each ticket's chance of winning is its weight, draw N winners with replacement, and let them be the new generation — all equal-weighted again.
The fat particle gets cloned three or four times; the starving particles die. You've concentrated your computation where the target actually lives. After the lottery you reset every weight to 1/N, so ESS jumps from near 1 back up to N.
But how you run the lottery matters. Multinomial resampling throws N independent random darts at a cumulative-weight ruler and keeps whichever particle each dart lands on. It works, but N independent darts inject a lot of Monte Carlo noise.
Systematic resampling is smarter: throw one random dart, then place the remaining N−1 darts at perfectly even spacing 1/N apart. One random offset, deterministic spacing. This injects strictly less variance than N independent draws while remaining unbiased — the expected number of offspring for each particle is still N·Wi.
The mechanism in both cases is the same: build the cumulative distribution of the weights (a staircase from 0 to 1), then for each pointer position u, find which step it falls in — that's the chosen parent's index. In code this is a single searchsorted into the CDF.
One counterintuitive fact, important enough to set the resampling schedule: resampling actually adds variance at the moment you do it. You're replacing distinct weighted particles with noisy clones, throwing away information. So you don't resample every step — you resample only when ESS drops below a threshold (typically N/2), and you prefer systematic resampling because it adds the least noise.
Top: weighted particles as unequal bars, with the cumulative-weight ruler (a staircase) beneath. Pointer marks drop onto the ruler; each mark selects the particle whose CDF bin it lands in. Bottom: the resampled swarm, now equal height, with clones stacked. Toggle Multinomial (N random marks) vs Systematic (one random offset, evenly spaced marks).
Normalized weights W = (0.1, 0.6, 0.2, 0.1). The cumulative CDF is the running sum: (0.1, 0.7, 0.9, 1.0). To remove randomness, we use given pointer positions.
Multinomial with fixed uniforms u = (0.05, 0.30, 0.50, 0.95). For each u, find the first CDF entry ≥ u:
Result: indices (0, 1, 1, 3). The fat particle 1 (weight 0.6) is chosen twice; particles 0 and 3 once each; particle 2 dies (multinomial happened to skip over its narrow 0.7–0.9 band).
Systematic with offset u0=0.1, N=4. The marks are (u0 + k)/N for k=0,1,2,3 = (0.1, 1.1, 2.1, 3.1)/4 = (0.025, 0.275, 0.525, 0.775). Searchsorted into the CDF:
Result: indices (0, 1, 1, 2) — the same fat particle is cloned, but the evenly-spaced marks pick up particle 2 that multinomial missed. After either scheme, all four weights reset to 1/4.
# systematic resampling: one offset, evenly spaced pointers import numpy as np def systematic_resample(W, rng): N = len(W) positions = (rng.random() + np.arange(N)) / N cum = np.cumsum(W / W.sum()) return np.searchsorted(cum, positions) # parent indices W = np.array([0.1, 0.6, 0.2, 0.1]) idx = systematic_resample(W, np.random.default_rng(0)) print(idx) # fat particle (idx 1) cloned
# library one-liner: then x = x[idx]; W = np.ones(N)/N idx = np.searchsorted(np.cumsum(W), (rng.random() + np.arange(N)) / N)
Your target is a nasty, sharply-peaked posterior — sampling it directly is hopeless. But a broad Gaussian is trivial to sample. What if you could deform the easy distribution into the hard one in a hundred tiny steps, carrying your particles along for the ride?
The geometric tempering path does exactly this. Take an easy start density p0 and the hard target pT, raise each to complementary powers, and multiply: πβ(x) ∝ p0(x)1−β · pT(x)β. As the inverse temperature β climbs from 0 to 1, the bridge slides from the pure easy world to the pure target.
Why this particular path? Because each intermediate distribution overlaps heavily with its neighbors. A particle never has to make a giant leap; one rung is only a little harder than the last. That overlap is what keeps the incremental weights well-behaved (next chapter) and lets a single cheap MCMC move (Chapter 8) re-equilibrate the swarm at each rung.
For two zero-mean Gaussians the geometry is gorgeous and exact. Multiplying Gaussian densities adds their log-densities, and a Gaussian's log-density is quadratic with coefficient equal to its precision (one over the variance). So the bridge precision is the linear interpolation of the two precisions: prec(β) = (1−β)·prec0 + β·precT.
This is the single most-confused point in tempering, so say it twice: the bridge interpolates precisions, not variances and not a mixture. A linear mixture of two Gaussians keeps two separate bumps and never narrows. The geometric path produces one Gaussian that smoothly contracts.
That smooth contraction is the whole point. Because no single rung is a cliff, particles flow continuously from the broad easy world into the sharp target, and the swarm stays populated the entire way down the staircase.
At β=0 the bridge is the broad easy Gaussian (std 3); at β=1 it is the sharp target (std 1); in between it contracts. Toggle a bimodal target to watch one broad blob split into two separated peaks as β rises (foreshadowing the showcase). The shaded band is the overlap between this rung and the next.
Easy p0 = N(0, 32), so prec0 = 1/9. Target pT = N(0, 12), so precT = 1. The bridge precision is prec(β) = (1−β)/9 + β.
The width shrinks smoothly from 3 to 1. Notice β=0.25 still has variance 3.0 — early rungs near the broad world barely move, which is why a good ladder spends more rungs where the action is. No single step is a cliff.
# geometric bridge for two zero-mean Gaussians: interpolate precision import numpy as np betas = np.linspace(0, 1, 5) for b in betas: prec = (1 - b) / 9 + b / 1 # 1/9 (easy) -> 1 (target) print(f'beta={b:.2f} var={1/prec:.3f} std={1/prec**0.5:.3f}')
# library one-liner: general unnormalized log-bridge log_bridge = lambda x, b: (1 - b) * logp0(x) + b * logpT(x) # geometric tempering
You've built the staircase of bridge distributions. Now the question: as a particle steps from rung βt−1 up to βt, how much should its weight change? The answer is beautiful and closed-form.
The weight update is the ratio of the new bridge to the old one, evaluated at the particle's position: winc = πβt(x) / πβt−1(x). For the geometric path this ratio collapses because the two bridges differ only in the exponents on the same two densities.
Write the bridge in log form: log πβ(x) = (1−β) log p0(x) + β log pT(x) (up to a constant). Subtract the two rungs and the p0 and pT terms regroup into a single difference scaled by the temperature step:
That's it — one exponential of (the temperature step) times (the log-density gap between target and base) at the particle's current x. Small β steps mean weights near 1, so the swarm barely degenerates per rung. This incremental weight is the atom from which both sampling and normalizing-constant estimation are built.
One critical subtlety: this weight is evaluated at the particle's position before the MCMC move that follows. The temperature step does the reweighting; the move does the exploring. In the standard SMC-sampler / AIS construction the move leaves the new bridge invariant and so adds nothing to the weight.
Conflating "reweight" and "move" is the classic SMC-sampler bug — it double-counts and breaks the unbiasedness of the evidence estimate. Keep them mentally separate: reweight = anneal (the temperature step), move = explore (the kernel).
The two faint curves are log p0 (broad base) and log pT (sharp target). Drag the particle across x: the canvas measures the vertical gap (log pT − log p0), multiplies by the slider's Δβ, and exponentiates to display the incremental weight as a bar. Large Δβ swings weights far from 1 (degeneracy risk); small Δβ keeps them near 1 (safe).
Unnormalized base f0(x) = exp(−x2/18) (the N(0,9) kernel) and target fT(x) = exp(−x2/2) (the N(0,1) kernel). Particle at x = 1.5, step Δβ = 0.25.
log f0(1.5) = −(1.52)/18 = −2.25/18 = −0.125.
log fT(1.5) = −(1.52)/2 = −2.25/2 = −1.125.
The gap: log fT − log f0 = −1.125 − (−0.125) = −1.000. The target likes x=1.5 less than the base does — it's a bit far out for std 1.
Scale by the temperature step: log winc = 0.25 × (−1.000) = −0.25. Exponentiate: winc = exp(−0.25) = 0.7788. The particle is gently down-weighted as the world sharpens around it.
# incremental weight for the geometric path import numpy as np def log_incremental_weight(x, b_prev, b, logp0, logpT): return (b - b_prev) * (logpT(x) - logp0(x)) # temp step * log-density gap logp0 = lambda x: -x**2 / 18 # N(0,9) kernel logpT = lambda x: -x**2 / 2 # N(0,1) kernel print(np.exp(log_incremental_weight(1.5, 0.0, 0.25, logp0, logpT))) # 0.7788
# library one-liner: add to running log-weight log_w_inc = (beta[t] - beta[t - 1]) * (logpT(x) - logp0(x))
Here is a thing standard MCMC can never tell you: the normalizing constant Z = ∫ f(x) dx, the total area under your unnormalized target. Metropolis-Hastings only ever needs ratios π(x′)/π(x), in which Z appears in both numerator and denominator and cancels exactly. The chain samples the shape but never learns the area underneath.
Yet Z is exactly the quantity you need most. In Bayesian model selection it is the model evidence p(data) — the marginal likelihood that tells you which model the data prefer. In energy-based models it is the partition function you need to turn an energy into a probability. No amount of mixing recovers it from a Metropolis chain.
SMC's incremental weights — the ones you were tempted to normalize away — secretly multiply up to the ratio ZT/Z0. The absolute, un-normalized weight mass is not garbage; averaged across the swarm and accumulated along the ladder, it reconstructs the ratio of normalizing constants between the target and the base.
The mechanism: at each rung the average of the incremental weights estimates the ratio πβt's normalizer to πβt−1's normalizer. The product of these per-rung ratios telescopes to ZT/Z0 — a ratio of areas, where MCMC could only ever see the shape.
And here's the payoff: if Z0 is a known easy distribution (a prior, a broad Gaussian) then you multiply by it and read off the absolute ZT for free. The base distribution is chosen precisely so its normalizer is known in closed form.
This is the unique selling point of the sampler view of SMC. Importance sampling, MCMC, and nested sampling each touch a piece of this, but SMC/AIS hand you samples from the hard target and its normalizing constant in one pass.
Easy p0 = N(0,9) with KNOWN area Z0 = √(2π·9); target pT = N(0,1) with area ZT = √(2π). As β steps up the ladder, the running product of mean incremental weights (the moving line) converges to the true ratio 1/3. The side bar contrasts: MCMC sees only the normalized shape; SMC tracks the absolute weight mass.
Base f0 = exp(−x2/18), an unnormalized N(0,9). Its true normalizer is Z0 = √(2π·9) = √(18π) = 7.5199.
Target fT = exp(−x2/2), an unnormalized N(0,1). Its normalizer is ZT = √(2π) = 2.5066.
The ratio: ZT/Z0 = 2.5066 / 7.5199 = 0.333333 = 1/3. This is no coincidence: for these zero-mean Gaussians Z scales with the standard deviation, so the ratio is just stdT/std0 = 1/3 = 0.333333.
This 1/3 is the number the SMC/AIS swarm reconstructs by accumulating incremental weights. And since Z0 = 7.5199 is known analytically, multiplying gives the absolute ZT = (1/3)·7.5199 = 2.5066 — the area MCMC could never compute.
# accumulate log incremental weights -> log(Z_T / Z_0) import numpy as np betas = np.linspace(0, 1, 200) logp0 = lambda x: -x**2 / 18; logpT = lambda x: -x**2 / 2 rng = np.random.default_rng(0); x = rng.normal(0, 3, 5000) # ~ p_0 logZ = 0.0 for t in range(1, len(betas)): lw = (betas[t] - betas[t - 1]) * (logpT(x) - logp0(x)) logZ += np.log(np.mean(np.exp(lw))) # (resample + MCMC move here in full SMC) print(np.exp(logZ)) # ~ Z_T / Z_0 = 1/3
# library one-liner: AIS / SMC evidence ratio logZ_ratio = sum(logsumexp(log_w_inc[t]) - np.log(N) for t in range(1, T))
Strip resampling out of SMC entirely and you get Radford Neal's Annealed Importance Sampling (AIS, 2001) — the workhorse behind every reported VAE log-likelihood and energy-based-model partition function in modern deep generative modeling.
Each particle is a lone climber. Start it in the easy distribution p0, ratchet β from 0 to 1, apply an MCMC move at each rung, and multiply in that rung's incremental weight (Chapter 5). No interaction between climbers, no resampling — just N independent ascents.
At the top, each climber hands you one importance weight for the target. Its total log-weight is simply the sum of its per-rung incremental log-weights — because log of a product is a sum of logs.
Now the crucial step. To estimate ZT/Z0 you average the weights in linear space, which in log space is the log-mean-exp: logsumexp(logw) − log(N). This is unbiased for log(ZT/Z0). AIS is embarrassingly parallel and is the standard tool for evaluating generative models.
The famous trap: do not take the mean of the log-weights. By Jensen's inequality, mean(log w) ≤ log(mean(w)), so averaging in log space systematically under-estimates the evidence. This is a real and famous bug when scoring VAEs — it makes a model look worse than it is.
AIS can still degenerate if the ladder is too coarse: the final-weight ESS warns you. The cure is the same dial as always — add rungs (finer ladder) until the incremental weights stay near 1 and the final ESS is healthy.
Several climbers ascend the β ladder simultaneously on separate tracks, each accumulating a running log-weight (the bar beside its track). No arrows between tracks — that's the no-resampling AIS signature. At the top, the per-particle log-weights pool into a log-mean-exp readout converging to log(ZT/Z0). Toggle resampling ON to morph back into full SMC.
First, one climber whose four per-rung log incremental weights are (−0.1, −0.2, −0.15, −0.05). Its total log-weight is the sum: −0.1 − 0.2 − 0.15 − 0.05 = −0.50.
Now pool four climbers with totals logw = (−0.5, −0.3, −0.8, −0.4). The AIS estimate is logsumexp(logw) − log(4).
Exponentiate each: e−0.5=0.6065, e−0.3=0.7408, e−0.8=0.4493, e−0.4=0.6703. Sum = 2.4669.
logsumexp = log(2.4669) = 0.9030. Subtract log(4) = 1.3863: estimate = 0.9030 − 1.3863 = −0.4833.
Compare to the naive mean of the log-weights, (−0.5−0.3−0.8−0.4)/4 = −0.50. The correct log-evidence estimate −0.4833 is higher — exactly the Jensen gap. The ESS of these four weights is 1/ΣWi2 = 3.8816 out of 4, a healthy run.
# AIS estimate of log(Z_T / Z_0) from per-particle log-weights import numpy as np from scipy.special import logsumexp def ais_logZ(particles_logw): N = len(particles_logw) return logsumexp(particles_logw) - np.log(N) # unbiased for log(Z_T/Z_0) logw = np.array([-0.5, -0.3, -0.8, -0.4]) print(ais_logZ(logw)) # -0.4833 (NOT the mean -0.5)
# library one-liner log_evidence_ratio = scipy.special.logsumexp(particle_logweights) - np.log(N)
Resampling clones the fat particle three times — but now you have three identical copies sitting on top of each other, which is useless for exploring. The fix: after resampling, take an MCMC step at the current temperature to jitter the clones apart.
You built the Metropolis acceptance rule in MCMC (that lesson owns detailed balance and why the chain leaves its target invariant). Here we only borrow it: the inter-rung kernel targets the current bridge πβ, so applying it re-spreads cloned particles without disturbing the rung's distribution.
The crucial twist is temperature softening. For a target with energy U(x) (so πβ ∝ exp(−βU(x))), the Metropolis acceptance of a proposed move is min(1, exp(−β·ΔU)), where ΔU is the energy change. The temperature β multiplies the barrier.
So warm rungs (small β) are forgiving: the same uphill move is accepted readily because β shrinks the effective barrier. Cold rungs (β near 1) are picky: the full barrier is felt and uphill moves are often rejected. This is precisely what lets particles cross between separated modes while the world is still warm.
And the same separation-of-concerns from Chapter 5 holds: this move changes the particle's position, not its weight. The kernel leaves πβ invariant, so it contributes nothing to the importance weight. Reweighting comes only from the temperature step.
Mixing these up — letting the move add to the weight — double-counts and breaks the unbiasedness of the Z estimate. Move = explore; reweight = anneal. Keep them separate and the machinery is exact.
A bimodal energy landscape with two wells. A particle proposes an uphill move; the acceptance gauge shows min(1, exp(−β·ΔU)). At WARM β you watch the particle hop the valley to the other mode (accepted because the bridge is flatter); at COLD β the same hop is usually rejected. Drag β and propose moves.
Energy U(x) = x2/2 (so the cold target is N(0,1)). A particle at x=1.0 proposes a move to x′=1.5 (uphill in energy).
U(1.0) = 1.02/2 = 0.5. U(1.5) = 1.52/2 = 1.125. So ΔU = U(x′) − U(x) = 1.125 − 0.5 = 0.625 (uphill).
Warm rung β=0.5: acceptance = min(1, exp(−0.5 × 0.625)) = exp(−0.3125) = 0.7316 — accept with probability 0.732.
Cold rung β=1.0: acceptance = min(1, exp(−1.0 × 0.625)) = exp(−0.625) = 0.5353 — accept with probability 0.535.
The same uphill move is far more likely to be accepted when the world is warm. That is how tempering lets particles roam between modes before the target sharpens and locks them in place.
# tempered Metropolis step: leaves pi_beta invariant, no weight change import numpy as np def metropolis_step(x, beta, U, step, rng): xp = x + step * rng.standard_normal() a = np.exp(-beta * (U(xp) - U(x))) # tempered acceptance return xp if rng.random() < a else x U = lambda x: x**2 / 2 print(np.exp(-0.5 * (U(1.5) - U(1.0)))) # 0.7316 warm print(np.exp(-1.0 * (U(1.5) - U(1.0)))) # 0.5353 cold
# library one-liner accept = np.exp(-beta * (U(x_prop) - U(x_cur))) # leave pi_beta invariant
Here's the payoff. A posterior with two well-separated modes — peaks at −4 and +4 with a deep valley between. Launch a single Metropolis chain at the right mode and it samples the right mode beautifully, reporting it as the whole answer, blind to the left mode for millions of steps because crossing the valley requires an astronomically rare uphill run.
Now launch the full SMC sampler: tempering ladder + resampling + MCMC moves, all three together. Start warm — one broad blob covering both peaks — then cool down the ladder. Particles distribute across both modes while the world is warm, resampling keeps both populated, and as β → 1 you land with the correct two-mode mixture.
The mechanism is everything you've built: Chapter 4's ladder deforms the broad base into the bimodal target; Chapter 5's incremental weights anneal each rung; Chapter 3's resampling fires when ESS drops below N/2; Chapter 8's temperature-softened moves let warm particles hop the valley before the barrier rises.
Why the single chain fails has a number behind it. Under the warm base N(0,9), the probability mass in the bridging region is large; under the cold target N(0,1) it is tiny. The warm world has a fat overlap connecting the modes; the cold world has a near-impassable gap.
Drive it yourself below. Widen the mode separation and watch the single chain fail harder while SMC keeps recovering both peaks. This is the showcase — no quiz; the simulation is the test.
The single broad base N(0,9) has std = √9 = 3.0, comparable to the half-separation 4 — particles spanning roughly [−3, 3] sit in the valley and can drift either way.
Compare the bridging-region mass. P(X>2) under the warm base N(0,9) is 1 − Φ(2/3) = 0.2525. The same P(X>2) under the cold target N(0,1) is 1 − Φ(2) = 0.0228.
The ratio: 0.2525 / 0.0228 = 11.08. There is over eleven times more mass in the bridging region when the world is warm — that fat overlap is exactly why warm particles cross the valley and cold ones get stranded.
Left: a lone Metropolis chain, stuck in one mode, histogram one-sided. Right: the full SMC swarm descending the β ladder (broad blob → two sharp peaks), with resample+move events and a histogram filling BOTH peaks near 50/50. The progress bar is β: 0→1; the ESS meter triggers resampling at N/2.
# full SMC sampler on a bimodal target (sketch) import numpy as np from scipy.special import logsumexp rng = np.random.default_rng(0) N, betas = 2000, np.linspace(0, 1, 50) logp0 = lambda x: -x**2 / 18 # broad base logpT = lambda x: logsumexp([-(x - 4)**2 / 2, -(x + 4)**2 / 2], axis=0) x = rng.normal(0, 3, N); logw = np.zeros(N) for t in range(1, len(betas)): lp = lambda b, x: (1 - b) * logp0(x) + b * logpT(x) logw += lp(betas[t], x) - lp(betas[t - 1], x) # incremental weight W = np.exp(logw - logsumexp(logw)) if 1 / np.sum(W**2) < N / 2: # resample on ESS x = x[rng.choice(N, N, p=W)]; logw[:] = 0 x = x + 0.5 * rng.standard_normal(N) # MCMC jitter (sketch) print((x > 0).mean()) # ~0.5: both modes found
# library one-liner: ladder + resample + kernel built in import blackjax; blackjax.tempered_smc(...) # or numpyro.infer SMC
You have two candidate models for the same dataset. Which one does the data prefer? The principled answer is the Bayes factor — the ratio of evidences p(data | model 1) / p(data | model 2). But each evidence is a brutal high-dimensional integral that MCMC can't touch.
AIS computes each one. Anneal from the prior (easy, with known normalizer) to the unnormalized posterior (the target), and the accumulated log-weights hand you log p(data) directly — the model evidence. This is literally how the field reports VAE log-likelihoods and EBM partition functions.
The pipeline below runs the whole thing: many climbers ascend a β ladder from prior to unnormalized posterior, each accumulating a log-weight; the live log-evidence (logsumexp − log N) stabilizes; a second model is annealed alongside; and a Bayes-factor needle swings to favor one model.
Watch the coarse-ladder bias live. With too few rungs, the incremental weights swing wildly, the log-space averaging bites (Jensen), and the evidence is under-estimated. Add rungs and the estimate converges to the true logZ overlaid for the tractable toy. The ESS meter warns you when the ladder is too coarse to trust.
This is the second showcase — no quiz; driving the pipeline to a correct, ESS-healthy evidence is the test. Read off the Bayes factor and render the verdict.
Two models scored by AIS return log-evidences logZ1 = −12.3 and logZ2 = −14.1 (each is log p(data | model k), the AIS output).
Log Bayes factor = logZ1 − logZ2 = −12.3 − (−14.1) = 1.8.
Bayes factor BF = exp(1.8) = 6.0496. The data are about 6× more probable under model 1 than model 2 — on the Jeffreys scale (BF 3–10 = "moderate"), this moderately favors model 1.
The pooling step under the hood: five climbers with logw = (−0.5, −0.3, −0.8, −0.4, −1.2) give logsumexp − log(5). Exponentials sum to 0.6065+0.7408+0.4493+0.6703+0.3012 = 2.7681; logsumexp = log(2.7681) = 1.0184; minus log(5) = 1.6094 gives −0.5912. Their ESS is 4.62 out of 5 — healthy, so trust the run. The whole point: each logZ came from AIS accumulating incremental weights, an integral MCMC could not compute.
Climbers ascend from prior to unnormalized posterior for two models. Live readouts: each model's log-evidence (logsumexp − log N), the true logZ overlay for the tractable toy, the ESS health check, and the log Bayes factor logZA − logZB with a verdict. Drop the rung count to watch the coarse-ladder bias under-estimate the evidence.
# AIS log-evidence and a Bayes factor for model selection import numpy as np from scipy.special import logsumexp def ais_log_evidence(log_prior_Z, particle_logw): # particle_logw = sum of incremental log-weights prior -> posterior return log_prior_Z + logsumexp(particle_logw) - np.log(len(particle_logw)) logZ1, logZ2 = -12.3, -14.1 # two models' AIS evidences print('log Bayes factor', logZ1 - logZ2) # 1.8 print('Bayes factor', np.exp(logZ1 - logZ2)) # 6.05 -> favor model 1
# library one-liner: VAE/EBM eval # anneal prior->posterior; logZ = logsumexp(logw) - log(N); BF = exp(logZ_A - logZ_B)
Step back and see the whole family. The particle filter moves a swarm through time to track a hidden state as data streams in — same particles, same resampling, but the sequence of distributions is the filtering posterior p(statet | data1:t), driven by physics.
The SMC sampler you just learned moves a swarm through a temperature ladder to reach one static target. The sequence is artificial, chosen by you, and the prize is samples from a hard distribution plus its normalizing constant. Same engine, opposite purpose. Knowing which one you need is half the battle.
SMC sits at the intersection of three lessons. It shares the importance weight w = π/q with importance sampling (Chapter 1). It shares the Metropolis/Gibbs kernel with MCMC (Chapter 8). It shares the particle-and-resample machinery with the particle filter (Chapter 3). And uniquely, it yields the normalizing constant.
For the gradient-powered kernels that make SMC moves efficient in high dimensions, the natural next step is Hamiltonian Monte Carlo — drop an HMC kernel into the inter-rung move and the swarm explores far faster. SMC and HMC compose cleanly.
The named historical references — Neal's AIS (2001), Del Moral, Doucet, and Jasra's SMC samplers (2006), Skilling's nested sampling — are still the backbone in 2024–2026. AIS is how modern papers report VAE and EBM log-likelihoods, and SMC ships in production probabilistic programming: NumPyro, Pyro, Stan's bridge-sampling, and BlackJAX's tempered SMC.
One diagnostic follows you everywhere: the weight ESS. It tells you whether to trust an evidence estimate, whether to add rungs, whether to resample. It is the through-line connecting every method in the Monte Carlo family.
| Method | Sequence of distributions | Gives Z? |
|---|---|---|
| Importance sampling | none (single proposal) | ratio only |
| MCMC | none (one fixed target) | no |
| Particle filter | data-driven state posterior over time | yes (likelihood) |
| SMC sampler | user-chosen tempering ladder, one static target | yes (evidence) |
| AIS | tempering ladder, no resampling | yes (evidence) |
Nodes for Importance Sampling, MCMC, Particle Filter, SMC Sampler, AIS, and HMC. Click a node to read its one-line relationship and the shared component. Toggle the split to color nodes by dynamic state (filtering) vs static target (sampling) — the boundary that makes the SMC-sampler-vs-particle-filter distinction vivid.
Given final AIS log-weights logw = (−0.5, −0.3, −0.8, −0.4, −1.2). Normalize: Wi = exp(logwi − logsumexp(logw)). From Chapter 10, logsumexp = 1.0184.
The normalized weights are exp(−0.5−1.0184)=0.2191, exp(−0.3−1.0184)=0.2676, exp(−0.8−1.0184)=0.1623, exp(−0.4−1.0184)=0.2421, exp(−1.2−1.0184)=0.1088 (summing to 1).
ESS = 1/ΣWi2 = 1/(0.0480+0.0716+0.0263+0.0586+0.0118) = 1/0.2164 = 4.62 out of 5. As a fraction: 4.62/5 = 0.924 — healthy, so trust this evidence estimate.
If ESS had collapsed toward 1, you'd add rungs before believing the logZ. This single diagnostic — the same ESS from Chapter 2 — is the through-line connecting every method in the family.
# the one diagnostic that follows you everywhere: weight ESS import numpy as np from scipy.special import logsumexp logw = np.array([-0.5, -0.3, -0.8, -0.4, -1.2]) W = np.exp(logw - logsumexp(logw)) print('ESS', 1 / np.sum(W**2), 'of', len(W)) # 4.62/5 -> trust the run
# next: hamiltonian-mc for gradient-powered kernels inside SMC moves # see: importance-sampling (weights), mcmc (kernels), particle-filter (filtering)