67. Importance sampling

Monte Carlo's error is \sigma/\sqrt n. The rate is fixed, so the only lever is \sigma — and importance sampling is the most powerful way to pull it.

It also solves a problem plain Monte Carlo simply cannot: estimating the probability of a rare event.

The failure it fixes

Estimate \theta = P(Z > 4) for a standard Normal. The true value is 3.167 \times 10^{-5}.

Draw 10,000 standard Normals and count how many exceed 4. The expected count is 10{,}000 \times 3.167\times 10^{-5} = 0.3.

Most runs produce zero hits, and the estimator confidently returns \hat\theta = 0.000000 — not "approximately zero", but exactly zero, with an estimated standard error of zero as well. The method reports certainty about a number it never sampled.

The problem is structural: almost every sample lands where h is zero, so almost every sample is wasted.

The identity

The fix is to sample from somewhere more useful and correct for the bias.

E_p[h(X)] = \int h(x)p(x)\,dx = \int h(x)\frac{p(x)}{q(x)}q(x)\,dx = E_q\left[h(X)\,w(X)\right]

with the importance weight

w(x) = \frac{p(x)}{q(x)}

The estimator becomes

\hat\theta = \frac{1}{n}\sum_{i=1}^n h(X_i)\,w(X_i), \qquad X_i \sim q

Sample from q, weight by p/q. The weights undo the distortion of having sampled from the wrong distribution — exactly, in expectation.

The requirement is a support condition: q(x) > 0 wherever h(x)p(x) \ne 0. If q never visits a region that matters, no weighting can recover it.

Contrast with rejection sampling

Both use a proposal q, and the difference is what they do with the mismatch:

• Rejection sampling discards proposals that don't fit, and the accepted ones are exact draws from p.
• Importance sampling reweights every proposal, keeping all of them.

Nothing is thrown away, which is why importance sampling can be dramatically more efficient — and why it needs no bounding constant M.

Choosing the proposal

The variance of the estimator is

\operatorname{Var}_q\big(h\,w\big) = E_q\big[h^2w^2\big] - \theta^2

Minimising it over q gives the optimal proposal

q^*(x) \propto |h(x)|\,p(x)

which achieves zero variance when h \ge 0.

That sounds too good, and it is: the normalising constant of q^* is \int h(x)p(x)dx = \theta — the very quantity you're trying to compute. The optimal proposal is unusable, but it tells you what to aim for:

Sample where |h| \cdot p is large — where the integrand actually contributes.

For a rare-event problem that means shifting the proposal into the tail.

The failure mode

Get this backwards and importance sampling is worse than useless.

If q has lighter tails than p, the weights w = p/q explode in the tail. A handful of samples end up with enormous weights, the estimator is dominated by them, and the variance can be infinite — while the sample standard error, computed from draws that never reached the tail, looks reassuringly small.

The rule: the proposal must have tails at least as heavy as the target. When in doubt, use something heavy-tailed like a t distribution.

A practical diagnostic is the effective sample size:

\text{ESS} = \frac{\left(\sum_i w_i\right)^2}{\sum_i w_i^2}

If 10,000 samples give an ESS of 12, your estimate rests on about 12 of them. Always check it.

Worked example

Estimate \theta = P(Z > 4) with 10,000 samples.

Plain Monte Carlo. Draw Z \sim N(0,1), count hits. Typically 0 hits, giving \hat\theta = 0 — a relative error of 100%.

Importance sampling. Take q = N(4, 1) — shifted so the region of interest is now the centre of the proposal. Roughly half the draws exceed 4.

The weight is a ratio of Normal densities:

w(x) = \frac{\varphi(x)}{\varphi(x - 4)} = \exp\left(-\frac{x^2}{2} + \frac{(x-4)^2}{2}\right) = \exp(8 - 4x)

Sample values and their weights:

x	w(x) = e^{8-4x}
4.0	3.355 \times 10^{-4}
4.5	4.540 \times 10^{-5}
5.0	6.144 \times 10^{-6}

Tiny weights — correctly so, because these are samples the target distribution would essentially never produce.

The result: relative error falls from 100% to around 2%, with a variance reduction of several thousand times, using the same 10,000 samples.

Note the tail condition holds: N(4,1) has the same Gaussian tail decay as N(0,1), so the weights e^{8-4x} shrink as x grows rather than exploding. Had we used a proposal with lighter tails, the same setup would have failed catastrophically.

Doing it in Python

The rare-event failure, and the fix:

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(0)
n, threshold = 10_000, 4.0
truth = 1 - norm.cdf(threshold)

# Plain Monte Carlo
plain = (rng.standard_normal(n) > threshold).mean()

# Importance sampling with q = N(4, 1)
X = rng.normal(threshold, 1, n)
w = norm.pdf(X) / norm.pdf(X, threshold, 1)
h = (X > threshold).astype(float)
imp = (h * w).mean()

print(f"true value        : {truth:.6e}")
print(f"plain Monte Carlo : {plain:.6e}")
print(f"importance sampling: {imp:.6e}")
print(f"\nplain relative error     : {abs(plain-truth)/truth:>8.2%}")
print(f"importance relative error: {abs(imp-truth)/truth:>8.2%}")

The variance reduction, measured across many runs:

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(1)
n, threshold, trials = 10_000, 4.0, 500
truth = 1 - norm.cdf(threshold)

plain, imp = [], []
for _ in range(trials):
    plain.append((rng.standard_normal(n) > threshold).mean())
    X = rng.normal(threshold, 1, n)
    w = norm.pdf(X) / norm.pdf(X, threshold, 1)
    imp.append(((X > threshold) * w).mean())

plain, imp = np.array(plain), np.array(imp)
print(f"{'method':>20} {'mean':>14} {'sd':>14} {'rel. error':>12}")
for name, e in [("plain", plain), ("importance", imp)]:
    print(f"{name:>20} {e.mean():>14.6e} {e.std():>14.6e} "
          f"{e.std()/truth:>11.2%}")

print(f"\nplain returned exactly 0 in {(plain == 0).mean():.1%} of runs")
print(f"variance reduction: {plain.var()/imp.var():,.0f}x")

Choosing the shift — the U-shaped variance curve:

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(2)
n, threshold, trials = 5_000, 4.0, 200
truth = 1 - norm.cdf(threshold)

print(f"{'proposal mean':>15} {'estimate':>14} {'sd of estimate':>16} {'ESS':>8}")
for mu in (0, 2, 3, 4, 5, 6, 8):
    ests, ess = [], []
    for _ in range(trials):
        X = rng.normal(mu, 1, n)
        w = norm.pdf(X) / norm.pdf(X, mu, 1)
        vals = (X > threshold) * w
        ests.append(vals.mean())
        ess.append(w.sum() ** 2 / (w ** 2).sum())
    ests = np.array(ests)
    print(f"{mu:>15} {ests.mean():>14.4e} {ests.std():>16.4e} {np.mean(ess):>8.0f}")

print(f"\ntrue value {truth:.4e}")
print("Best near mu = 4 -- the region that matters. Too small wastes samples;")
print("too large makes the weights erratic.")

The catastrophic failure — a proposal with lighter tails:

import numpy as np
from scipy.stats import norm, t as t_dist

rng = np.random.default_rng(3)
n, trials = 20_000, 300

# Target: a t distribution with 3 df (HEAVY tails). Estimate E[|X|].
target = t_dist(3)
truth = 2 * target.expect(lambda x: x, lb=0)

for name, sampler, density in [
    ("Normal proposal (LIGHT tails)", lambda k: rng.normal(0, 2, k),
     lambda x: norm.pdf(x, 0, 2)),
    ("t(2) proposal (heavy tails)", lambda k: rng.standard_t(2, k),
     lambda x: t_dist(2).pdf(x)),
]:
    ests, ess = [], []
    for _ in range(trials):
        X = sampler(n)
        w = target.pdf(X) / density(X)
        ests.append((np.abs(X) * w).mean())
        ess.append(w.sum() ** 2 / (w ** 2).sum())
    ests = np.array(ests)
    print(f"{name:>32}: mean {ests.mean():>7.4f}  sd {ests.std():>8.4f}  "
          f"ESS {np.mean(ess):>7.0f}")

print(f"\ntrue E[|X|] = {truth:.4f}")
print("\nThe light-tailed proposal gives erratic estimates and a collapsed ESS:")
print("its weights explode in the tail, where the target still has mass.")

Your turn

1. Why does plain Monte Carlo fail for rare events?

2. What is the optimal proposal, and why can't you use it?

3. Your importance weights range from 10^{-8} to 10^{6}. What's wrong?

Solutions

1. Because almost every sample contributes zero.

Estimating \theta = P(X \in A) for a rare A means averaging an indicator that is nearly always 0. With n samples the expected number of hits is n\theta, so for \theta = 10^{-5} and n = 10^4 you expect 0.3 hits.

Two consequences:

• The estimate is often exactly 0, which is not merely imprecise but qualitatively wrong — it reports impossibility.
• The relative error is terrible. For an indicator, \operatorname{Var} = \theta(1-\theta) \approx \theta, so

\frac{\operatorname{SE}}{\theta} \approx \frac{1}{\sqrt{n\theta}}

To reach 10% relative error you need n\theta \approx 100, i.e. $n \approx 100/\theta$. For \theta = 10^{-6} that's 100 million samples.

Importance sampling changes the problem so that most samples land in A and carry information.

2. The variance-minimising proposal is

q^*(x) \propto |h(x)|\,p(x)

and for h \ge 0 it gives exactly zero variance — every weighted sample returns the same value.

You can't use it because its normalising constant is

\int |h(x)|p(x)\,dx

which for h \ge 0 is precisely \theta, the integral you set out to compute. Knowing q^* requires already knowing the answer.

It's still useful as a target to approximate. It says: put proposal mass where |h| \cdot p is large. In a rare-event problem that's the tail; in a Bayesian problem it's the region of high posterior mass. Adaptive importance sampling and cross-entropy methods iteratively fit a tractable q toward q^*.

3. The proposal has tails that are too light relative to the target, and the estimate is not trustworthy.

Weights spanning 14 orders of magnitude mean a few samples carry essentially all the mass:

\text{ESS} = \frac{(\sum w_i)^2}{\sum w_i^2}

will be tiny — perhaps 1 or 2 — so 20,000 draws are effectively a couple of draws. And in the worst case \operatorname{Var}(hw) is infinite, in which case the CLT doesn't apply and the reported standard error is meaningless.

The insidious part: the sample standard error looks fine, because the samples that would reveal the problem are the ones q almost never produces. The diagnostic must be the ESS, not the standard error.

Fixes:

• Use a heavier-tailed proposal — a t distribution is the standard choice.
• Defensive mixture: q_{\text{new}} = 0.9q + 0.1p, which bounds the weights by 10.
• Weight clipping / truncation, which trades a little bias for a large variance reduction.
• Check the ESS every time, before believing any importance-sampling result.

Check yourself in code

Show that plain Monte Carlo fails on a rare event while importance sampling succeeds, using the same number of samples.

Print exactly this:

truth 3.167e-05
plain 0.000e+00
importance 3.105e-05
importance is closer: True

Use default_rng(0), n = 10000, threshold 4.0, and a proposal of N(4, 1). Format every estimate in scientific notation with 3 decimal places.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(0)
n, threshold = 10_000, 4.0
truth = 1 - norm.cdf(threshold)

print(f"truth {truth:.3e}")

plain = (rng.standard_normal(n) > threshold).mean()
print(f"plain {plain:.3e}")

# Sample from N(4,1), weight by norm.pdf(X)/norm.pdf(X, 4, 1), and report the
# importance-sampling estimate plus whether it is closer to the truth.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(0)
n, threshold = 10_000, 4.0
truth = 1 - norm.cdf(threshold)

print(f"truth {truth:.3e}")

plain = (rng.standard_normal(n) > threshold).mean()
print(f"plain {plain:.3e}")

X = rng.normal(threshold, 1, n)
w = norm.pdf(X) / norm.pdf(X, threshold, 1)
imp = ((X > threshold) * w).mean()
print(f"importance {imp:.3e}")
print("importance is closer:", bool(abs(imp - truth) < abs(plain - truth)))

---

Importance sampling draws from a proposal q and corrects with weights w = p/q, so no sample is wasted. Aim the proposal where |h|p is large — the zero-variance optimum q^* \propto |h|p is unusable, since its normaliser is the answer, but it tells you the direction. Keep the proposal's tails at least as heavy as the target's, and check the effective sample size, because a proposal that is too light produces exploding weights and a standard error that lies.

---

That completes the course. You began with nine salaries and the question of what "typical" means, and you've arrived at the computational methods that make modern statistics possible — by way of the axioms, random variables, the limit theorems that justify inference, estimation and testing, regression, both schools of inference, stochastic processes, the measure-theoretic foundations, and information theory.

The thread running through all of it: probability describes uncertainty, and statistics reasons backwards from data to the process that produced it. Everything else is machinery for doing that honestly.