Importance Sampling: reweight
the right samples

Chapter 1: Reweight — Sample Here, Answer There

Your problem in Chapter 0 was that p almost never visits the dangerous region. So let's cheat — deliberately. We draw samples from a proposal q that we get to choose, one that does visit the tail, by pushing its mean toward the danger zone.

But now the samples are biased toward danger: they over-represent the tail and under-represent the bulk. If we just averaged them we would get the wrong answer. The fix is pure bookkeeping. Each sample x carries an honest correction factor — the importance weight w(x) = p(x)/q(x) — that says "q over-produced this point, so down-weight it."

The magic is that this correction is exact, not approximate. Watch the algebra. We want an expectation under p, and we multiply and divide by q inside the integral:

Read it left to right: the expectation of f under p equals the expectation of w(x)f(x) under q. So we can sample from q, compute w(x)f(x) for each draw, average, and get E_p[f]. This is the importance sampling identity, and it is the entire foundation of the lesson.

The only requirement is that q is nonzero wherever p·f is nonzero (you can't reweight a region you never sample). Subject to that, the identity holds for any valid q. The estimator is unbiased no matter where q sits. What q changes is not the mean — it's the variance, which is the subject of the next chapter.

Play with the widget: drag q's mean and width. Dots are drawn from q, so they cluster under q's hump. Each dot's size encodes its weight w = p/q — dots in q's over-sampled region shrink (down-weighted), dots in p-favored regions grow. The weighted average tracks the true E_p[f] no matter where you drag q.

Let's compute weights by hand to demystify them. Take p = N(0,1) and q = N(3,1). Both are Gaussians, so the ratio of their densities collapses to a clean exponential — the normalizing constants cancel and only the exponents remain:

Now evaluate at three fixed points (no randomness — this is deterministic arithmetic). At x = 2.5: exponent = −3(2.5) + 4.5 = −7.5 + 4.5 = −3, so w = e⁻³ = 0.049787. At x = 3.5: −10.5 + 4.5 = −6, so w = e⁻⁶ = 0.002479. At x = 4.0: −12 + 4.5 = −7.5, so w = e^−7.5 = 0.000553.

Every weight is far below 1, and they shrink fast as x grows. That makes sense: q is centered at 3, so it over-represents the deep tail relative to p (centered at 0), and the correction pulls those over-sampled points sharply down. The weight is doing exactly its job — undoing q's bias point by point.

Here is the whole estimator in code. We sample around 3 yet recover E_p[X] ≈ 0.

python — importance sampling from scratch
import numpy as np
from scipy.stats import norm
rng = np.random.default_rng(0)
N = 100_000
x = rng.normal(3.0, 1.0, N)               # draw from proposal q
w = norm.pdf(x, 0, 1) / norm.pdf(x, 3, 1)   # p/q correction
f = x                                      # estimating E_p[X] (should be ~0)
est = np.mean(w * f)
print(est)                                   # ~0 despite sampling around 3

python — the whole estimator is one line
est = np.mean(w * f)   # a weighted average, nothing more

Chapter 2: Variance — The Two-Faced Estimator

You just learned IS is unbiased for any q. So why not pick q sloppily? Because unbiased does not mean accurate. Run the same IS estimator with two different proposals and one converges in a hundred samples while the other still flails after a million. The mean is right either way; the spread is what kills you.

Importance sampling is a sharp knife. With a proposal that hugs where the integrand has its mass, it cuts variance dramatically below crude MC. With a proposal that misses that mass, it cuts you — the variance can be larger than crude MC, even literally infinite. Same estimator, two faces.

The variance of the IS estimator is the variance of the random quantity w(x)f(x) under q, divided by N. So we care about the second moment E_q[(w f)²]. Writing it out exposes the danger:

Look at where q sits: in the denominator. Wherever q is tiny but p·f is not, the integrand blows up. If q has lighter tails than p (it decays faster), then far out in the tail p/q → ∞, this integral can diverge, and the IS variance is infinite even though the estimate stays unbiased. That is the whole story of IS variance in one line.

Run the widget to feel it. One slider morphs q from "matched" to "mismatched." Watch the right-hand running estimate destabilize: it sits flat, then leaps when one giant weight finally appears. The weight histogram develops a fat tail with a single dominating spike — the visual signature of exploding IS variance.

Let's anchor the safe case with exact arithmetic. Estimate E_p[f] with p = N(0,1) and f(x) = x². The truth is E_p[x²] = 1 (the variance of a standard normal). Take the "do-nothing" proposal q = p, so w = 1 everywhere and IS reduces to plain MC.

The per-sample variance of plain MC here is Var(x²) = E[x⁴] − (E[x²])². The fourth moment of a standard normal is E[x⁴] = 3 (a fact you can look up or derive from the moment generating function). So:

That 2 is the baseline. A heavier-tailed q (larger width) keeps w = p/q bounded and can push the per-sample variance below 2. A lighter-tailed q makes w blow up exactly where x² is largest, so E_q[(w x²)²] can diverge — the variance goes to infinity while the mean stays glued to 1.

The practical takeaway: in IS you watch the variance, not the mean. A stable-looking running mean that suddenly jumps is your alarm bell.

python — same mean, wildly different variance
import numpy as np
from scipy.stats import norm
rng = np.random.default_rng(1)
N = 200_000; f = lambda x: x**2
def is_est(q_sigma):
    x = rng.normal(0, q_sigma, N)
    w = norm.pdf(x,0,1)/norm.pdf(x,0,q_sigma)
    vals = w*f(x)
    return vals.mean(), vals.var()   # mean ~1 always; var explodes if q_sigma<1
print(is_est(1.5))   # safe (heavier tails)
print(is_est(0.5))   # dangerous (lighter tails)

python — watch var, not mean
mean, var = (w*f(x)).mean(), (w*f(x)).var()

Chapter 3: Optimal q — The Proposal That Beats Noise

If choosing q changes variance, what is the best possible q? Intuitively, put your samples exactly where the integrand f(x)p(x) is large — don't waste a single draw on regions that contribute nothing. Push that intuition to its limit and something startling falls out.

The variance-minimizing proposal for estimating E_p[f] (with f ≥ 0) is

That is, q* is shaped like the product |f|·p, renormalized to integrate to 1. Now watch what happens to the weighted integrand under q*. We substitute q* into w(x)f(x) = (p/q*)·f and the p and f factors cancel:

The weighted integrand is a constant — it equals Z for every x. A constant random variable has zero variance. So under q*, a single sample returns the exact answer Z with no error at all. This is the zero-variance proposal, and it is the north star that tells us which way to push q.

Make it concrete with p = N(0,1) and f(x) = x². Then Z = E_p[x²] = 1, and q*(x) = x²·φ(x)/1. Check the cancellation at two arbitrary points. At x = 2: q*(2) = 4·φ(2), so w = φ(2)/(4φ(2)) = 1/4, and w·f = (1/4)·4 = 1. At x = 0.5: q*(0.5) = 0.25·φ(0.5), so w = 1/0.25 = 4, and w·f = 4·0.25 = 1. Both give exactly Z = 1.

Notice how the weight and the integrand trade off perfectly: where f is large the weight is small, where f is small the weight is large, and the product is pinned at Z everywhere. That balance is the entire point.

Play with the widget: a blend slider morphs q from p toward q* ∝ |f|p. As you approach q*, the per-sample w·f bars flatten to a single level and the variance readout collapses toward zero. Then flip the "reveal the catch" toggle.

So the practical doctrine is: you cannot have q*, but you can imitate it. Match the shape of |f|p, and crucially keep tails at least as heavy as p·|f| (Chapter 2's warning) so the imitation never blows up.

python — the integrand is constant under q*
import numpy as np
from scipy.stats import norm
# f = x^2, Z = E_p[x^2] = 1, so q* = x^2 * phi(x) / Z
Z = 1.0
for x in [0.5, 2.0, -1.3]:
    p = norm.pdf(x); qstar = (x**2 * p) / Z
    w = p/qstar; print(x, w*x**2)   # all == Z exactly

python — the recipe is guidance, not code
# pick q with shape ~ |f(x)|*p(x), and tails >= p*|f|

Chapter 4: ESS — How Many Samples Do You Really Have?

You drew 1000 weighted samples — but one of them has weight 100 and the rest have weight 1. Is that really 1000 samples' worth of information? Obviously not: the weighted average is dominated by that single sample. You would be a fool to trust it like an honest 1000.

There is a one-line diagnostic that converts "N weighted samples" into "how many fair, equal-weight samples this is actually worth." It is the effective sample size, and it is the single most important health check in all of importance sampling.

Read the formula's behavior at its two extremes. If all N weights are equal, the numerator is (Nw)² = N²w² and the denominator is Nw², so ESS = N — every sample counts fully. If one weight dwarfs all others, the sums are both dominated by that one term, the ratio collapses toward 1 — you effectively have a single sample.

Why this exact formula? It comes from setting the variance of a weighted mean of equal-information samples equal to the variance you actually have. The ratio (∑w)²/∑w² is precisely the "equivalent number of equal samples" — it is the inverse of the sum of squared normalized weights, which is a standard concentration measure.

Let's compute it by hand for two fixed weight vectors of length N = 5 — pure arithmetic, no randomness. First a moderately skewed vector w = [0.1, 0.2, 0.05, 4.0, 0.3].

Sum: 0.1 + 0.2 + 0.05 + 4.0 + 0.3 = 4.65, so (∑w)² = 4.65² = 21.6225. Sum of squares: 0.01 + 0.04 + 0.0025 + 16.0 + 0.09 = 16.1425. Therefore ESS = 21.6225 / 16.1425 = 1.3395. Five skewed samples are worth barely more than one — that single weight of 4.0 eats the set.

Now a degenerate vector w = [1, 1, 1, 1, 100]. Sum: 4 + 100 = 104, so (∑w)² = 10816. Sum of squares: 4×1 + 100² = 4 + 10000 = 10004. So ESS = 10816 / 10004 = 1.0812 — essentially one effective sample despite having five. The lone giant weight swallows everything.

Drag the skew slider in the widget and watch the ESS dial swing from N (healthy, even weights) down toward 1 (degenerate, one spike). The efficiency readout ESS/N tells you what fraction of your sampling effort is actually paying off.

In practice you compute ESS every run and treat a low value as a stop sign — re-center q, widen it, or reach for a more powerful method. We'll use ESS again in Chapter 6 (over-tilting), Chapter 9 (too few annealing rungs), and Chapter 10 (the universal diagnostic).

python — effective sample size from scratch
import numpy as np
def ess(w):
    w = np.asarray(w, float)
    return w.sum()**2 / (w**2).sum()
print(ess([0.1,0.2,0.05,4.0,0.3]))   # 1.3395
print(ess([1,1,1,1,100]))            # 1.0812 -> degenerate

python — equivalent one-liner
ess = (w.sum()**2) / (w**2).sum()   # or 1/np.sum((w/w.sum())**2)

Chapter 5: Self-norm — When You Only Know p Up to a Constant

In Bayesian inference your posterior is p(θ|data) = likelihood × prior / evidence, and the evidence in the denominator is a brutal integral you usually cannot compute. So you can write down the posterior only up to that mystery constant. You know the shape, not the scale.

Standard IS needs the true normalized p(x) inside the weight w = p/q. If you only have an unnormalized p̃(x) = c·p(x) for some unknown c, ordinary IS gives you c·E_p[f] — off by the very factor you can't compute. Dead end? No. There is a trick so clean it feels like a cancellation by magic.

Form the weights from the unnormalized density, then divide the weighted sum by the sum of the weights. The unknown constant appears in both the numerator and the denominator, so it cancels. This is self-normalized importance sampling (SNIS).

Because w̃ = c·w, both sums scale by c, and c divides itself out. As a free bonus, the denominator mean(w̃) is itself an estimate of the unknown normalizer — we recover Ẑ without ever computing the evidence integral directly.

Let's do it by hand with four fixed samples. Take unnormalized target p̃(x) = exp(−x²/2)·exp(x), which is N(1,1) times an unknown constant. Use proposal q = N(0,1). The normal kernels cancel cleanly: w̃(x) = p̃(x)/q(x) = exp(−x²/2 + x)/exp(−x²/2) = exp(x).

Evaluate the unnormalized weights at the fixed samples x = [−1, 0, 1, 2]: w̃ = [e⁻¹, e⁰, e¹, e²] = [0.3679, 1.0, 2.7183, 7.3891]. Their sum is 0.3679 + 1.0 + 2.7183 + 7.3891 = 11.4753.

Now the weighted sum of x: (−1)(0.3679) + (0)(1.0) + (1)(2.7183) + (2)(7.3891) = −0.3679 + 0 + 2.7183 + 14.7782 = 17.1286. Divide:

The true mean of N(1,1) is 1; with only four samples we land at 1.49, in the right ballpark — and crucially we never needed the normalizer of p̃, because it cancelled in the ratio. With more samples this converges to 1 exactly.

In the widget, a slider scales p̃ by an arbitrary constant c (mimicking the unknown normalizer). Ordinary IS shifts with c — broken. SNIS stays pinned on the truth — visual proof the constant cancels. The recovered Ẑ tracks c on its own.

This single trick is what makes IS viable across Bayesian posteriors, statistical physics (where the partition function is the unknown constant), and variational inference. Almost every real target you meet is unnormalized; SNIS is how you cope.

python — self-normalized IS from scratch
import numpy as np
from scipy.stats import norm
rng = np.random.default_rng(2)
x = rng.standard_normal(50_000)            # proposal q = N(0,1)
log_ptilde = -x**2/2 + x                    # unnormalized log target (N(1,1)*const)
w = np.exp(log_ptilde - norm.logpdf(x))     # w-tilde = exp(x)
est = np.sum(w * x) / np.sum(w)             # SNIS -> ~1
print(est)

python — a weighted mean IS self-normalized IS
est = np.average(f(x), weights=w)   # numpy does the ratio for you

Chapter 6: Rare events — Tilting Toward Disaster

Back to the catastrophe estimate from Chapter 0 — the one that returned a useless 0. Now you have the machinery: the IS identity (Chapter 1), the variance warning (Chapter 2), the q* north star (Chapter 3), and ESS to keep you honest (Chapter 4). It is time to actually estimate a rare-event probability.

Don't sample from the loss model and pray. Sample from a tilted model — a proposal shifted so that most samples land in the danger zone — then reweight each by w = p/q to undo the tilt. A million-to-one event now shows up in every batch, and your estimate finally has an error bar.

This is the killer application of importance sampling. In structural reliability you tilt toward the failure surface; in finance you tilt toward the Value-at-Risk threshold; in telecom you tilt toward buffer overflow. The recipe is always the same: shift mass into the rare region, then correct.

Let's estimate P(X > 3) for X ∼ N(0,1), whose true value is 0.0013499. Tilt the proposal to q = N(3,1), centered exactly on the threshold so that roughly half the samples cross it. The weight is the same closed form we derived in Chapter 1: w(x) = exp(−3x + 4.5).

Use four fixed proposal samples x = [2.5, 3.2, 3.5, 4.0] and the indicator f = 1{x > 3}. First the weights. w(2.5) = exp(−7.5 + 4.5) = exp(−3) = 0.049787. w(3.2) = exp(−9.6 + 4.5) = exp(−5.1) = 0.006097. w(3.5) = exp(−10.5 + 4.5) = exp(−6) = 0.002479. w(4.0) = exp(−12 + 4.5) = exp(−7.5) = 0.000553.

Now the indicators. Only x = 2.5 fails the test x > 3, so 1{x > 3} = [0, 1, 1, 1]. Multiply weights by indicators: w·f = [0, 0.006097, 0.002479, 0.000553]. The IS estimate is the average:

From just four samples we land at 0.00228, the same order of magnitude as the truth 0.00135 — whereas crude MC would need around 7400 samples just to see ten hits. The tilt-then-correct trick has bought us a factor of roughly two thousand in efficiency, and more for rarer events.

In the widget, slide q's mean from 0 toward and past the threshold. Watch the IS estimate's error bar shrink dramatically as q approaches the threshold — while crude MC stays stuck at 0/N. Push q too far past t and the ESS readout drops, teaching the Goldilocks lesson below.

This connects directly to Chapter 3: the optimal proposal for a tail probability concentrates right at the boundary of the rare set, where the indicator switches on but p still has appreciable density. Centering q at the threshold is a cheap, excellent imitation of q*.

python — rare-event probability by tilting
import numpy as np
from scipy.stats import norm
rng = np.random.default_rng(3)
N = 10_000; t = 3.0
x = rng.normal(3.0, 1.0, N)              # tilt q toward the tail
w = norm.pdf(x,0,1)/norm.pdf(x,3,1)      # = exp(-3x+4.5)
est = np.mean(w * (x > t))               # ~0.00135 with a real error bar
print(est, norm.sf(t))

python — with the closed-form weight inlined
est = np.mean(np.exp(-3*x+4.5) * (x > 3))   # tilted estimate of P(X>3)

Chapter 7: Off-policy — Grading a Policy You Never Ran

A hospital has months of logged treatment decisions made under its current protocol — the behavior policy b. You have designed a new target policy π and you want its expected outcome without running it on real patients. You cannot re-collect data under π; that experiment is exactly what you're trying to avoid.

This is off-policy evaluation, and importance sampling is the classical answer. A trajectory is just a sequence of states and actions; the probability of a trajectory factors over its steps. Reweighting from b to π means multiplying by the ratio of trajectory probabilities — and because the environment dynamics cancel (they're shared), only the policy ratios remain.

You met the per-step ideas in policy gradients and Q-learning; here the new idea is the product of per-step ratios and why it is a variance time bomb.

Here G is the trajectory's observed return and ρ is the cumulative importance ratio. Multiply the return by ρ and you get an unbiased estimate of the target policy's value from a behavior-policy trajectory. The danger lives entirely in that product symbol.

Let's compute it for a 3-step trajectory. Along the actions actually taken, the target probs are π = [0.8, 0.5, 0.9] and the behavior probs are b = [0.4, 0.5, 0.3], with observed return G = 2.0. The per-step ratios are π/b = [0.8/0.4, 0.5/0.5, 0.9/0.3] = [2.0, 1.0, 3.0].

The cumulative weight is the product: ρ = 2.0 × 1.0 × 3.0 = 6.0. The IS return estimate is ρ·G = 6.0 × 2.0 = 12.0. So this one logged trajectory, which earned a return of 2, counts as evidence for a target-policy value of 12 — because π was six times more likely to take this exact path than b was.

Now the variance bomb, with a tiny toy. Suppose each step's ratio equals 2 or 0 with equal probability (mean 1, so still unbiased). Then E[ratio²] per step = 0.5·4 + 0.5·0 = 2. Over H = 10 independent steps, E[ρ²] = 2¹⁰ = 1024, so Var(ρ) = 1024 − 1 = 1023.

The second moment grows exponentially in the horizon H. Double the trajectory length and the variance roughly squares. This is the curse of long-horizon off-policy evaluation, and it is why importance sampling alone struggles on anything but short episodes.

In the widget, a horizon slider lengthens the trajectory and a mismatch slider widens the gap between π and b. Watch the spread of cumulative weights across many sampled trajectories balloon as either grows. Toggle to weighted IS (WIS) to see its tighter but slightly biased spread.

The latest research lives here. Modern off-policy evaluation for LLM agents, recommender systems, and clinical policies leans on doubly-robust estimators (2024–2026) that combine an IS correction with a learned value model: the value model handles the bulk, IS corrects the residual, and the variance is a fraction of pure IS. Weighted IS and per-decision IS are the building blocks underneath them.

python — off-policy IS return estimate
import numpy as np
# logged trajectory: action probs under behavior b and target pi, plus return G
pi = np.array([0.8, 0.5, 0.9]); b = np.array([0.4, 0.5, 0.3]); G = 2.0
rho = np.prod(pi / b)            # cumulative importance ratio = 6.0
is_return = rho * G             # ordinary IS estimate = 12.0
print(rho, is_return)

python — weighted IS across many trajectories
value_hat = np.average(returns, weights=rhos)   # self-normalized, lower variance

Chapter 8: Renderer — Sampling the Light That Matters

You're rendering a glossy floor. Each pixel's color is an integral over all incoming light directions, weighted by how the surface reflects them. The reflectance function is the BRDF (bidirectional reflectance distribution function), and the pixel value is the integral of incoming radiance times the BRDF times a cosine term over the hemisphere.

Sample directions uniformly over the hemisphere and most rays fly off where the surface barely reflects anything — wasted samples, and the image is a blizzard of noise. Importance sampling says: draw directions from the reflection lobe itself, so every ray lands where light actually contributes. This is THE reason every offline film frame is renderable at all.

This is importance sampling in disguise. The "proposal" q is the pdf you sample directions from; the "weight" is the estimator f(dir)/pdf(dir). Matching pdf to the integrand is exactly Chapter 3's q* idea, now in a graphics costume.

Consider a Lambertian (matte) surface with constant incoming radiance L = 1. The pixel integral is the integral over the hemisphere of L·cos(θ)/π. The single-sample estimator is f(dir)/pdf(dir). Let's compare two pdfs at a direction with cos(θ) = 0.8.

The integrand there is f = L·cos/π = 0.8/π = 0.2546. Under uniform sampling the pdf is 1/(2π) = 0.15915, so the estimator is f/pdf = 0.2546/0.15915 = 1.6. At a grazing direction cos = 0.2 the same uniform estimator is (0.2/π)/(1/(2π)) = 0.4. Uniform estimates swing from 0.4 to 1.6 across directions — high variance.

Under cosine-weighted sampling the pdf is cos/π = 0.8/π = 0.25465, so the estimator is f/pdf = 0.2546/0.25465 = 1.0. At cos = 0.2 it is (0.2/π)/(0.2/π) = 1.0 again — in fact it is 1.0 at every direction, because the cos in the integrand cancels the cos in the pdf.

A constant estimator has zero variance: for constant radiance, cosine-weighted sampling renders the exact pixel value from a single sample. This is the textbook win of BRDF importance sampling — the proposal perfectly cancels the integrand, just like q* in Chapter 3.

This is the showcase. Toggle uniform vs cosine sampling; sharpen the lobe with the glossiness slider (making good sampling matter more); sweep the sample count. Watch the two pixel patches: the importance-sampled one converges smoothly to clean while the uniform one stays grainy at the same budget. The f/pdf bars tell the whole story — uniform's vary wildly, cosine's are flat.

python — uniform vs cosine estimator variance
import numpy as np
rng = np.random.default_rng(4)
N = 20_000; L = 1.0
# uniform hemisphere: pdf = 1/(2pi)
c_u = rng.random(N)                       # cos(theta) proxy ~ uniform
est_u = (L*c_u/np.pi)/(1/(2*np.pi))    # f/pdf, high variance
# cosine-weighted: pdf = cos/pi -> estimator is L, constant
c_c = np.sqrt(rng.random(N))             # cosine-weighted cos(theta)
est_c = (L*c_c/np.pi)/(c_c/np.pi)         # all == L
print(est_u.var(), est_c.var())          # cosine variance ~ 0

python — the one line every path tracer runs
# wi = brdf.sample(); weight = brdf(wi)*cos / pdf(wi)   # f/pdf at the sampled direction

No quiz here — the simulation is the test. Toggle the sampler back and forth at low sample count and watch which patch you would actually ship.

Chapter 9: Annealed IS — Bridging Two Worlds With a Chain of Proposals

Your proposal q and target p are so different that a single p/q reweighting has near-zero ESS — every weight is either 0 or astronomical. From Chapter 4 you know that means you effectively have one sample no matter how many you draw. More draws don't help when q and p barely overlap.

The fix is a relay race. Don't jump from q to p in one leap; lay down a ladder of intermediate distributions that morph gradually from the easy q to the hard p, and let importance weights accumulate one rung at a time. This is annealed importance sampling (AIS), and it is the doorway to the entire Sequential Monte Carlo lesson next door.

The intermediate distributions are usually geometric blends: p_β ∝ p₀^(1−β)·p₁^β as β sweeps from 0 to 1. At each rung you nudge a particle with an MCMC move targeting p_β, then multiply in an incremental weight — the ratio of consecutive unnormalized densities. The accumulated weight estimates the normalizer ratio Z₁/Z₀.

Let's see the telescoping with a single deterministic particle. Take base p₀ = N(0,1) with normalizer Z₀ = 1, and unnormalized target p̃₁(x) = exp(−x²/2)·exp(x), which is N(1,1) up to a constant. The true normalizer ratio is the Gaussian integral:

Now hold a particle at x = 0.5 (no MCMC move, for illustration) and let the per-rung density ratios chain. The intermediate-density ratios telescope — consecutive numerators and denominators cancel — so the product collapses to p̃₁(x)/p₀(x) = exp(x) regardless of how many rungs you used.

At x = 0.5 the accumulated weight is exp(0.5) = 1.6487, which is exactly the true normalizer ratio Z₁/Z₀ = exp(0.5) = 1.6487. So the average of accumulated AIS weights recovers the otherwise-intractable normalizer ratio — the "evidence" or partition-function ratio that Bayesian model comparison and IWAE-style bounds both crave.

This is the showcase. The widget shows a ladder of densities morphing from p₀ to p₁; a particle walks down it; the accumulated log-weight grows rung by rung; and the running Z₁/Z₀ estimate and ESS update live. Crank the "rungs" slider from 1 (single-step IS) upward and watch ESS climb as the estimate stabilizes onto exp(0.5).

This is a bridge, not the destination. The full unbiased log-Z estimator — averaging accumulated weights with the log-mean-exp machinery — belongs to Sequential Monte Carlo, which generalizes AIS by adding resampling between rungs. Here we just plant the seed: chain IS reweightings with MCMC moves along a tempering ladder.

python — annealed importance sampling (sketch)
import numpy as np
rng = np.random.default_rng(5)
def logp0(x): return -x**2/2
def logp1(x): return -x**2/2 + x          # unnormalized target
betas = np.linspace(0,1,8)                  # the annealing ladder
x = rng.standard_normal(10_000); logw = np.zeros_like(x)
for b0,b1 in zip(betas[:-1], betas[1:]):
    logw += (b1-b0)*(logp1(x)-logp0(x))     # incremental AIS weight
    # (here: + an MCMC move targeting p_b1 between rungs)
Zratio = np.exp(logw).mean()                # -> ~exp(0.5) = 1.6487
print(Zratio)

python — the bridge to SMC
# AIS = a chain of IS reweightings with MCMC moves; the SMC lesson adds resampling

No quiz here — the simulation is the test. Set rungs to 1 with a wide gap and watch ESS collapse; then add rungs and watch it recover.

Chapter 10: Connections — Where the Weights Lead

You started Chapter 0 with a useless 0 / 1,000,000. You now hold a precision instrument: sample where it's convenient, reweight to where it's true, and always check ESS to know whether you can trust the answer. That single discipline — w = p/q, then watch ESS — is the whole lesson.

But every chapter also hinted at a wall. Long horizons make the off-policy weight product explode (Chapter 7). Near-disjoint q and p collapse the ESS to one (Chapter 9). Lighter-tailed proposals send the variance to infinity (Chapter 2). Those walls are exactly where importance sampling hands off to its more powerful cousins.

Here is the map. Importance sampling is the reweighting primitive sitting underneath rare-event simulation, off-policy RL, modern variational inference, and sequential Monte Carlo. Hover the nodes in the widget to surface each one's idea and its characteristic failure mode; the same ESS diagnostic routes you between them.

Let's close with a capstone consistency check that ties three chapters together with one shared number. Take the degenerate weight vector from Chapter 4, w = [1, 1, 1, 1, 100], and ask what it implies for every downstream application — rare-event, off-policy, AIS. They all read the same diagnostic.

We already computed ESS = (∑w)²/∑w² = 104²/10004 = 1.0812. The efficiency is ESS/N = 1.0812/5 = 0.2162, i.e. about 22% — these five weighted samples carry the information of barely one.

Whether this weight vector came from a rare-event tilt that overshot (Chapter 6), an off-policy trajectory product (Chapter 7), or a too-coarse AIS ladder (Chapter 9), the verdict is identical: the proposal or bridge is degenerate, trust nothing, and add structure. Re-center q, switch to weighted IS or a doubly-robust estimator, or add AIS rungs and move to SMC. One diagnostic, every application.

The threads forward: Sequential Monte Carlo generalizes the annealing ladder of Chapter 9; policy gradients and Q-learning use the off-policy ratios of Chapter 7 inside their updates; and Monte Carlo is the root where it all began.

python — the one diagnostic every IS-based method shares
import numpy as np
def ess(w):
    w = np.asarray(w, float); return w.sum()**2 / (w**2).sum()
# rare-event tilt, off-policy products, AIS ladders -> all read ESS
w = np.array([1,1,1,1,100.0])
print(ess(w), ess(w)/len(w))   # 1.08, 22% efficiency -> degenerate

python — the universal stop sign
if ess(w)/len(w) < 0.1: print('proposal degenerate -> add structure (SMC/WIS/DR)')