Monte Carlo: rejection sampling, importance sampling and variance reduction

Before this lesson: lesson 02 (expectation and inverse-CDF sampling), lesson 03 (the exponential and normal distributions), lesson 06 (the law of large numbers, the CLT and the Monte Carlo error bar), lesson 10 (posteriors known only up to a constant). From foundations: log-sum-exp. Time: about 75 minutes, plus the notebook.

---

Why this matters #

• Expectations everywhere. A Bayesian prediction is $\mathbb{E}_{p(\theta \mid \mathcal D)}[p(y \mid \theta)]$; a VAE trains on a Monte Carlo estimate of the ELBO; a policy gradient is a Monte Carlo average over trajectories.
• Importance weights in modern ML. PPO's probability ratio, off-policy evaluation in reinforcement learning, the IWAE bound and particle filters are all importance sampling. Speculative decoding in LLM serving is a form of rejection sampling.
• Rare events in science and engineering. The probability that a flood, an earthquake or a structural load exceeds a design level is tiny. Plain sampling needs about $100/\mu$ draws for 10% accuracy on an event of probability $\mu$; importance sampling can cut this by orders of magnitude.
• The bridge to MCMC. When no good proposal exists (high dimension, complex posterior), Markov chain Monte Carlo takes over (lesson 15). Its error bars use the same ideas as here.

1. Predict first #

2. The picture #

Estimate $\pi$ with darts. Throw $n$ points uniformly into the square $[-1, 1]^2$ and count the fraction inside the unit disk. The disk covers $\pi/4$ of the square, so $4 \times$ (fraction inside) estimates $\pi$. Each dart is a coin flip with $P(\text{inside}) = \pi/4$, so the estimate has standard error $4\sqrt{(\pi/4)(1 - \pi/4)/n} \approx 1.64/\sqrt n$.

Interactive figure: loading…

The error formula has three levers. The standard error of an average of $n$ i.i.d. draws of $f(X)$ is $\sigma_f/\sqrt n$, where $\sigma_f$ is the standard deviation of $f(X)$. To make it smaller you can:

• take more samples. This is slow: 4 times the work buys 2 times the accuracy.
• sample from somewhere else. Draw from a proposal $q$ that puts points where they matter, and correct with weights $w = p/q$ (importance sampling).
• make the integrand less variable. Subtract a correlated quantity with a known mean (control variates), or pair each draw with a mirrored twin (antithetic variables).

Weights correct the bias. Suppose the target is $p = \mathcal N(0, 1)$ but we draw from a wider $q$. Points far out are over-represented, so they get a small weight $p(x)/q(x)$; points near 0 are under-represented and get a weight above 1. The weighted average is then correct on average. But if $q$ is narrower than $p$, a rare draw in the tail gets an enormous weight and dominates the whole average.

Interactive figure: loading…

3. Definitions and main results #

Simple Monte Carlo (recap of lesson 06) #

We want $\mu = \mathbb{E}_p[f(X)] = \int f(x)\,p(x)\,dx$. Draw $X_1, \dots, X_n$ i.i.d. from $p$ and set

Then $\mathbb{E}[\hat\mu_n] = \mu$ and $\operatorname{Var}(\hat\mu_n) = \sigma_f^2/n$ with $\sigma_f^2 = \operatorname{Var}_p(f(X))$. The law of large numbers gives $\hat\mu_n \to \mu$, and the CLT gives the error bar $\hat\mu_n \pm 1.96\,s_f/\sqrt n$, where $s_f$ is the sample standard deviation of the $f(X_i)$ (lesson 06).

Rare events. For a probability $\mu = P(A)$, take $f = \mathbf 1_A$. Then $\sigma_f^2 = \mu(1 - \mu)$ and the relative standard error is

Here $n\mu$ is the expected number of hits. For 10% relative error you need about 100 hits, so about $100/\mu$ draws. The error depends on how many times you see the event, not on $n$ alone.

Exact samplers: inverse CDF and rejection #

Importance sampling is easier to understand after two methods that produce exact draws.

Inverse CDF (lesson 02). If $U \sim \text{Uniform}(0, 1)$ and $F$ is a CDF, then $F^{-1}(U)$ has CDF $F$ (with the generalised inverse $F^{-1}(u) = \inf\{x : F(x) \ge u\}$ for steps and flat parts). For the exponential with rate $\lambda$, $x = -\log(1 - u)/\lambda$. It needs a CDF you can invert; the Gaussian has none in closed form (use Box–Muller, lesson 05, or a numerical inverse).

Rejection sampling. Often you can evaluate the target only up to a constant: $p = \tilde p/Z$, where you know $\tilde p$ but not $Z = \int\tilde p$. A posterior is the typical case: likelihood times prior is $\tilde p$, and $Z$ is the evidence (lesson 10). Choose a proposal density $q$ that you can sample, and a constant $M$ with an envelope: $\tilde p(x) \le M q(x)$ for all $x$.

Algorithm: draw $X \sim q$ and $U \sim \text{Uniform}(0, 1)$ independently. Accept $X$ if $U \le \tilde p(X)/(Mq(X))$; otherwise try again.

Picture: the points $(X, U\,Mq(X))$ are uniform under the graph of $Mq$. We keep the ones under the graph of $\tilde p$. Those are uniform under $\tilde p$, so their $x$-coordinate has density proportional to $\tilde p$.

Theorem (rejection sampling is exact). An accepted $X$ has density $p$. Each trial is accepted with probability $Z/M$, so the number of trials per accepted sample is geometric with mean $M/Z$.

Proof.

1. Given $X = x$, the acceptance probability is $\tilde p(x)/(Mq(x))$, a number in $[0, 1]$ because of the envelope.
2. For a set $A$: $P(X \in A,\ \text{accept}) = \int_A q(x)\,\frac{\tilde p(x)}{Mq(x)}\,dx = \frac1M\int_A \tilde p(x)\,dx$.
3. Take $A$ = everything: $P(\text{accept}) = Z/M$.
4. Divide: $P(X \in A \mid \text{accept}) = \frac{(1/M)\int_A\tilde p}{Z/M} = \int_A p(x)\,dx$. So the accepted $X$ has density $p$.
5. Trials are i.i.d. with success probability $Z/M$, so their number is geometric with mean $M/Z$. $\blacksquare$

Use the smallest valid $M = \sup_x \tilde p(x)/q(x)$. The supremum must be finite, so $q$ must have tails at least as heavy as $\tilde p$: a Gaussian proposal can never envelope a Cauchy target. If $M$ is too small, the envelope dips below $\tilde p$ and the samples are silently wrong there. A too large $M$ is only slow.

Interactive figure: loading…

High dimension kills rejection. For the target $\mathcal N(0, I_d)$ and the proposal $\mathcal N(0, s^2 I_d)$ with $s > 1$, the ratio $p/q$ is largest at $x = 0$, so $M = s^d$ and the acceptance rate is $s^{-d}$. With $s = 1.1$ and $d = 100$ it is $7 \cdot 10^{-5}$. Rejection sampling is a tool for low dimension.

Importance sampling #

Importance sampling keeps every draw from $q$ and fixes the bias with weights.

The identity. If $q(x) > 0$ wherever $f(x)p(x) \ne 0$, then with the importance weight $w(x) = p(x)/q(x)$,

Estimator. Draw $X_1, \dots, X_n$ i.i.d. from $q$ and set $\hat\mu_q = \frac1n\sum_i w(X_i) f(X_i)$. It is unbiased by the identity. It is an average of i.i.d. terms, so

Estimate the standard error by the sample standard deviation of the terms $w(X_i)f(X_i)$, divided by $\sqrt n$.

The best proposal. $q^* \propto |f|\,p$ minimises $\sigma_q^2$.

1. Only $\mathbb{E}_q[w^2 f^2]$ depends on $q$.
2. By Jensen (or Cauchy–Schwarz), $\mathbb{E}_q[w^2 f^2] \ge \big(\mathbb{E}_q[w|f|]\big)^2 = \big(\mathbb{E}_p|f|\big)^2$ for every $q$.
3. Equality holds when $w|f|$ is constant, that is when $q \propto |f|\,p$. $\blacksquare$

For $f \ge 0$ the best proposal $q^* = fp/\mu$ gives zero variance, but it needs the unknown $\mu$. The practical lesson: put the proposal where $|f|\,p$ is large. For the rare event $\{X > a\}$, that region is just above $a$, so shift the proposal there.

Interactive figure: loading…

Why the proposal needs heavier tails. Take $p = \mathcal N(0, 1)$ and $q = \mathcal N(0, s^2)$. Then $w(x) = s\exp\!\big(-\tfrac{x^2}{2}(1 - s^{-2})\big)$ and

When $s < 1$, the weight grows without bound in the tails: $q$ rarely visits them, and a visit gets a huge weight. A heavier-tailed $q$ ($s > 1$) has bounded weights, $w \le s$. In general the condition is $\int f^2 p^2/q < \infty$.

Infinite variance is a silent failure. The estimator stays unbiased, and it still converges (the strong law of large numbers needs only $\mathbb{E}_q|wf| < \infty$). But there is no CLT error bar, and the sample standard error is itself unreliable. Most runs miss the tail and look stable; for $f \ge 0$ they are usually too small. A rare run jumps.

Self-normalised importance sampling #

Often you know only $\tilde p = Zp$. Use unnormalised weights $w_i = \tilde p(x_i)/q(x_i)$ and

This is self-normalised importance sampling (SNIS).

• $\hat Z$ is unbiased for $Z$, because $\mathbb{E}_q[\tilde p/q] = \int\tilde p = Z$. For a Bayesian model this is the evidence $p(\mathcal D)$. Its log is biased low: $\mathbb{E}[\log\hat Z] \le \log Z$ by Jensen.
• $\tilde\mu$ is consistent: the numerator over $n$ tends to $Z\mu$ and the denominator over $n$ to $Z$.
• $\tilde\mu$ is biased, because $\mathbb{E}[A/B] \ne \mathbb{E}[A]/\mathbb{E}[B]$. With $n = 1$, $\tilde\mu = f(X_1)$ with $X_1 \sim q$, so $\mathbb{E}[\tilde\mu] = \mathbb{E}_q[f]$, not $\mathbb{E}_p[f]$. The bias is $O(1/n)$, small next to the standard error $O(1/\sqrt n)$.
• Its variance is approximately $\frac1n\mathbb{E}_q\big[w^2(f - \mu)^2\big]$ (normalised $w = p/q$; delta method). Estimate it with $\sum_i\bar w_i^2\big(f(x_i) - \tilde\mu\big)^2$.
• Numerics. Work with log-weights $\ell_i = \log\tilde p(x_i) - \log q(x_i)$ and compute $\bar w_i = \exp\big(\ell_i - \operatorname{LSE}(\ell)\big)$ with $\operatorname{LSE}(\ell) = m + \log\sum_i e^{\ell_i - m}$, $m = \max_i\ell_i$ (log-sum-exp). Raw weights overflow or underflow.

The effective sample size #

Definition (Kish). $\displaystyle n_e = \frac{\big(\sum_i w_i\big)^2}{\sum_i w_i^2} = \frac{1}{\sum_i \bar w_i^2} \in [1, n].$

Where it comes from. Take i.i.d. $Y_i$ with variance $\sigma^2$ and fixed weights. The weighted mean $\sum_i w_iY_i/\sum_i w_i$ has variance $\sigma^2\sum_i w_i^2/(\sum_i w_i)^2$. An unweighted mean of $n_e$ draws has variance $\sigma^2/n_e$. Set the two equal and solve for $n_e$.

Properties. $n_e = n$ for equal weights, and $n_e \to 1$ when one weight dominates. Scaling all weights by a constant does not change it, so $\tilde p$ is enough. Its population version is $n/\mathbb{E}_q[w^2]$ (normalised $w$). For $p = \mathcal N(0, 1)$ and $q = \mathcal N(0, 2^2)$ this is $n\sqrt 7/4 \approx 0.661\,n$.

The curse of dimension. For $p = \mathcal N(0, I_d)$ and $q = \mathcal N(m, I_d)$ one finds $\mathbb{E}_q[w^2] = e^{\|m\|^2}$, so the population ESS is $n\,e^{-\|m\|^2}$. A shift of 0.1 standard deviations in each of $d = 1000$ coordinates gives $\|m\|^2 = 10$ and an ESS of $4.5 \cdot 10^{-5}\,n$. Small mismatches multiply over coordinates, and the weights collapse onto one draw.

A diagnostic, not a guarantee. The ESS uses only the weights, not $f$. A small ESS proves a problem. A large ESS does not prove success: if $q$ misses a mode of $p$ entirely, all draws can have similar weights and the estimate is still wrong.

Variance reduction #

Control variates. Let $g$ be a function with a known mean $\gamma = \mathbb{E}[g(X)]$. For any constant $\beta$,

is unbiased, because the subtracted term has mean 0.

1. Its variance is $\frac1n\big(\sigma_f^2 - 2\beta\operatorname{Cov}(f, g) + \beta^2\sigma_g^2\big)$, a parabola in $\beta$.
2. Set the derivative to zero: $\beta^* = \operatorname{Cov}(f, g)/\operatorname{Var}(g)$, the least-squares slope of $f$ on $g$.
3. Substitute: the minimum is $(1 - \rho^2)\,\sigma_f^2/n$, where $\rho = \operatorname{Corr}(f, g)$. $\blacksquare$

With $\rho = 0.9$ the variance falls to 19%; the sign of $\rho$ does not matter. In practice $\beta$ is estimated from the same draws, which adds a bias of order $1/n$. It helps when $g$ is cheap, strongly correlated with $f$, and has an exactly known mean. This is the same idea as CUPED in A/B testing, where a pre-experiment metric is the control (lesson 12), and as the baseline in a policy gradient.

Antithetic variables. Pair each uniform $U$ with $1 - U$ (or a Gaussian $\varepsilon$ with $-\varepsilon$). With $n$ evaluations of $h$, use $n/2$ pairs and average $Y_j = \frac12\big(h(U_j) + h(1 - U_j)\big)$.

1. $\operatorname{Var}(Y_j) = \frac14\big(\sigma^2 + \sigma^2 + 2\rho\sigma^2\big) = \frac{\sigma^2}{2}(1 + \rho)$, with $\rho = \operatorname{Corr}\big(h(U), h(1 - U)\big)$.
2. So the estimator has variance $\frac{\sigma^2(1 + \rho)/2}{n/2} = \sigma^2(1 + \rho)/n$, against $\sigma^2/n$ for $n$ independent draws.
3. It helps exactly when $\rho < 0$. $\blacksquare$

If $h$ is monotone, $h(U)$ and $h(1 - U)$ move in opposite directions, so $\rho \le 0$ and antithetics never hurt. If $h(u) = h(1 - u)$ (symmetric about $\frac12$), then $\rho = 1$ and the variance doubles. If $h$ is linear, the variance is 0.

Interactive figure: loading…

Answers to section 1. (1) About 50%. The expected number of hits is $n\mu \approx 3.2$, so the relative standard error is $\sqrt{(1 - \mu)/(n\mu)} \approx 0.56$. With probability $e^{-3.2} \approx 0.04$ there is no hit at all. (2) Proposal A. Its variance is $s^2 = 0.25 < 1/2$, so $\mathbb{E}_q[w^2] = \infty$: the estimate is unbiased but has infinite variance. Proposal B has bounded weights $w \le 2$ and an ESS of about $0.661\,n$.

4. Worked example #

(a) Rejection: a half-normal from an exponential. Target $\tilde p(x) = e^{-x^2/2}$ on $[0, \infty)$, with $Z = \sqrt{\pi/2} = 1.2533$. Proposal $q(x) = e^{-x}$ (Exponential(1)).

• $\tilde p/q = e^{x - x^2/2}$. Its exponent $x - x^2/2$ has derivative $1 - x$, so the maximum is at $x = 1$: $M = e^{1/2} = 1.6487$.
• Accept $x$ with probability $\tilde p(x)/(Mq(x)) = e^{-x^2/2 + x - 1/2} = e^{-(x - 1)^2/2}$.
• Acceptance rate $Z/M = 1.2533/1.6487 = 0.760$; on average $1/0.760 = 1.32$ proposals per sample.
• Proposal $x = 0.5$ with $u = 0.6$: $e^{-0.125} = 0.882 > 0.6$, accept. Proposal $x = 2.5$ with $u = 0.6$: $e^{-1.125} = 0.325 < 0.6$, reject.

(b) Importance sampling by hand. Target $p = \mathcal N(0, 1)$, proposal $q = \mathcal N(0, 2^2)$, $f(x) = x^2$ (true value 1). Then $w(x) = p(x)/q(x) = \frac{e^{-x^2/2}/\sqrt{2\pi}}{e^{-x^2/8}/(2\sqrt{2\pi})} = 2e^{-3x^2/8}$. Four draws from $q$:

$x_i$	$-2$	$0$	$1$	$3$
$w_i$	0.446	2.000	1.375	0.068
$w_i f(x_i)$	1.785	0	1.375	0.616

• Plain IS: $\frac14(1.785 + 0 + 1.375 + 0.616) = 0.944$.
• SNIS: $\sum w_i = 3.889$, so $\tilde\mu = 3.776/3.889 = 0.971$; and $\hat Z = 3.889/4 = 0.972$ (true $Z = 1$).
• ESS: $\sum w_i^2 = 6.093$, so $n_e = 3.889^2/6.093 = 2.48$ of 4 draws. The draw at 0 alone carries 51% of the normalised weight.

(c) A rare event. $\mu = P(X > 4) = 3.167 \cdot 10^{-5}$ for $X \sim p = \mathcal N(0, 1)$, $f(x) = \mathbf 1\{x > 4\}$. Proposal $q = \mathcal N(4, 1)$.

• The weight: $\log p(x) - \log q(x) = -\frac{x^2}{2} + \frac{(x - 4)^2}{2} = -4x + 8$, so $w(x) = e^{-4x + 8}$. On the event it is small: $w(4) = e^{-8} = 3.4 \cdot 10^{-4}$.
• Second moment: $\mathbb{E}_q[w^2 f] = \mathbb{E}_p[wf] = \int_4^\infty e^{-4x + 8}\varphi(x)\,dx$, where $\varphi$ is the standard normal density. Complete the square: $e^{-4x}\varphi(x) = e^{8}\varphi(x + 4)$, so the integral is $e^{16}P(X > 8) = 5.53 \cdot 10^{-9}$.
• Variance per draw: IS $5.53 \cdot 10^{-9} - \mu^2 = 4.52 \cdot 10^{-9}$; plain MC $\mu(1 - \mu) = 3.17 \cdot 10^{-5}$.
• Ratio about 7,000. With $n = 10^5$, the relative standard error falls from 0.56 to 0.0067.

(d) Variance reduction for $\ln 2 = \mathbb{E}\big[\frac{1}{1 + U}\big]$. Let $h(u) = 1/(1 + u)$.

• $\mathbb{E}[h^2] = \int_0^1 (1 + u)^{-2}du = \frac12$, so $\operatorname{Var}(h(U)) = \frac12 - (\ln 2)^2 = 0.01955$.
• Control $g(U) = U$ (mean $\frac12$, variance $\frac1{12}$). $\mathbb{E}[Uh(U)] = \int_0^1\frac{u}{1 + u}du = 1 - \ln 2$, so $\operatorname{Cov}(h, U) = (1 - \ln 2) - \frac{\ln 2}{2} = -0.0397$. Then $\beta^* = -0.0397 \cdot 12 = -0.477$ and $\rho^2 = \frac{0.0397^2}{0.01955/12} = 0.969$. The variance falls by the factor 0.031: about 32 times fewer draws.
• Antithetic: $\mathbb{E}[h(U)h(1 - U)] = \int_0^1\frac{du}{(1 + u)(2 - u)} = \frac23\ln 2$, so $\operatorname{Cov} = \frac23\ln 2 - (\ln 2)^2 = -0.0184$ and $\rho = -0.939$. The factor is $1 + \rho = 0.061$: about 16 times fewer evaluations. $h$ is monotone, so $\rho < 0$ was guaranteed.

5. In code #

The notebook asks you to write the samplers and estimators yourself. Here NumPy and SciPy check the numbers of sections 3 and 4 by simulation, with $n = 10^5$ draws each.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 100_000

# Simple Monte Carlo: pi with a CLT error bar
u = rng.uniform(-1, 1, size=(n, 2))
f = 4 * (np.sum(u**2, axis=1) <= 1)
print(f"{f.mean():.4f} +- {2 * f.std(ddof=1) / np.sqrt(n):.4f}")            # 3.1471 +- 0.0104

# Rejection: half-normal target exp(-x^2/2) on [0, inf), Exp(1) proposal, M = e^{1/2}
x = rng.exponential(size=n)
keep = rng.uniform(size=n) <= np.exp(-(x - 1)**2 / 2)                     # p~ / (M q)
print(round(keep.mean(), 3), round(stats.kstest(x[keep], stats.halfnorm.cdf).pvalue, 2))   # 0.757 0.62

# Importance sampling for P(X > 4) with q = N(4, 1)
mu = stats.norm.sf(4)
plain = (rng.standard_normal(n) > 4).astype(float)
xq = 4 + rng.standard_normal(n)
wf = np.exp(-4 * xq + 8) * (xq > 4)                                       # w = p/q times the indicator
for name, v in [("plain", plain), ("IS", wf)]:
    print(name, f"{v.mean():.3e} rel. SE {v.std(ddof=1) / np.sqrt(n) / mu:.4f}")
print(f"exact {mu:.3e}")
# plain 5.000e-05 rel. SE 0.7060      (5 hits)
# IS 3.176e-05 rel. SE 0.0067
# exact 3.167e-05

# Self-normalised IS with log-weights: E[X^2] = 1 under p = N(0, 1), q = N(0, 2^2)
xs = 2 * rng.standard_normal(n)
logw = stats.norm.logpdf(xs) - stats.norm.logpdf(xs, scale=2)
wbar = np.exp(logw - logw.max()); wbar /= wbar.sum()
print(round(np.sum(wbar * xs**2), 3), round(1 / np.sum(wbar**2) / n, 3))   # 0.997 0.661

# Control variate and antithetic pairs for ln 2 = E[1/(1+U)], same number of evaluations
U = rng.uniform(size=n)
h = 1 / (1 + U)
beta = np.cov(h, U)[0, 1] / np.var(U, ddof=1)
cv = h - beta * (U - 0.5)
anti = 0.5 * (1 / (1 + U[: n // 2]) + 1 / (2 - U[: n // 2]))             # n/2 pairs = n evaluations
print(round(beta, 3), round(cv.var() / h.var(), 3), round(2 * anti.var() / h.var(), 3))   # -0.477 0.031 0.061

In this run plain Monte Carlo saw 5 hits instead of the expected 3.2, so its estimate is 58% too high; its sample standard error is itself noisy. Importance sampling is within 0.3% of the exact value.

6. Where it shows up #

• Speculative decoding in LLMs (ml-systems 12). A small draft model proposes a token $x \sim q$; the large model accepts it with probability $\min(1, p(x)/q(x))$. On rejection it samples once from the normalised residual $\max(0, p - q)$. The accepted part has mass $\min(p, q)$ and the residual part $\max(0, p - q)$; they add up to $p$, so the output distribution is exactly the large model's.
• Importance weights in reinforcement learning. Off-policy evaluation reweights returns by $\pi(a \mid s)/\mu(a \mid s)$ (RL 04). PPO's probability ratio is an importance weight that is clipped to control its variance (RL 09).
• Variational inference. The IWAE bound $\mathbb{E}\big[\log\frac1K\sum_k w_k\big]$ is $\mathbb{E}[\log\hat Z]$ with $K$ draws from the encoder: biased low by Jensen, and tighter as $K$ grows (generative models 01).
• Particle filters and data assimilation. A particle filter is sequential SNIS: weight the ensemble members by the likelihood of new observations, and resample when the ESS drops (a common rule is $n_e < n/2$). In high-dimensional weather models the weights collapse onto one member unless the ensemble grows exponentially with the effective dimension. This is one reason weather centres use ensemble Kalman and variational methods instead (ml-for-science 06).
• Rare events in science. Seismic hazard (the probability that ground motion at a site exceeds a design level), extreme rainfall, failure of a structure. A proposal shifted toward the event, as in (c), cuts the cost by orders of magnitude.
• Control variates. The baseline in REINFORCE (RL 07, calculus 12). Multi-fidelity Monte Carlo uses a cheap emulator, whose mean is known from many cheap runs, as the control for an expensive simulator.
• Antithetic sampling. Mirrored noise $(\varepsilon, -\varepsilon)$ in evolution strategies and in finite-difference gradient estimators.

7. Common mistakes #

• A proposal with lighter tails than the target. The weights have infinite variance; the estimate looks stable for many runs and then jumps. Check $\int f^2p^2/q < \infty$, or use a heavy-tailed proposal (a Student $t$, a wider Gaussian, or a mixture with a broad component).
• Forgetting the support condition. If $q = 0$ where $fp \ne 0$, plain IS is biased. SNIS needs $q > 0$ wherever $p > 0$.
• Weights in linear space. $\tilde p(x)/q(x)$ overflows or underflows in high dimension. Use log-weights and log-sum-exp.
• Trusting a large ESS. It cannot see a mode that no draw visited.
• Using the plain-IS standard error for SNIS. Use $\sum_i\bar w_i^2(f(x_i) - \tilde\mu)^2$, not the sample variance of $w_if(x_i)$.
• A rejection envelope that is too low. If $\tilde p > Mq$ somewhere, the samples are silently wrong there. A too high $M$ is only slow.
• Reporting a rare-event estimate with too few hits. With 3 expected hits, the estimate is often 0 or off by a factor 2. Count the hits; aim for at least 100, or use importance sampling.
• Antithetic pairs for a symmetric integrand. For $h(u) = h(1 - u)$ they double the variance.
• A control variate whose mean is only estimated. The error of that estimate adds directly to the error of $\hat\mu_\beta$; $\mathbb{E}[g]$ must be known exactly, or far more precisely than $\mu$.

8. Check yourself #

1. You sample $\mathcal N(0, 1)$ by rejection from a standard Cauchy proposal $q(x) = 1/(\pi(1 + x^2))$. Find $M$ and the acceptance rate. Can you sample a Cauchy target from a Gaussian proposal?
   Answer

   $p/q = \sqrt{\pi/2}\,(1 + x^2)e^{-x^2/2}$. The derivative of $(1 + x^2)e^{-x^2/2}$ is $x(1 - x^2)e^{-x^2/2}$, so the maximum is at $x = \pm 1$: $M = \sqrt{\pi/2}\cdot 2e^{-1/2} = \sqrt{2\pi}\,e^{-1/2} \approx 1.520$. Acceptance rate $1/M \approx 0.658$ (normalised target, $Z = 1$). The reverse is impossible: Cauchy/Gaussian $\to \infty$ in the tails, so no finite $M$ exists.
2. How many plain Monte Carlo draws do you need to estimate $P(X > 3) = 1.35 \cdot 10^{-3}$ for $X \sim \mathcal N(0, 1)$ with 10% relative standard error?
   Answer

   Relative SE $\approx 1/\sqrt{n\mu} = 0.1$ gives $n\mu = 100$, so $n \approx 100/(1.35 \cdot 10^{-3}) \approx 74{,}000$ draws (the exact formula $(1 - \mu)/(0.01\mu)$ gives about 73,980).
3. Show that SNIS is biased with a one-line example. Why is it still used?
   Answer

   With $n = 1$, $\tilde\mu = f(X_1)$ with $X_1 \sim q$, so $\mathbb{E}[\tilde\mu] = \mathbb{E}_q[f] \ne \mathbb{E}_p[f]$ in general. It is used because it needs $p$ only up to a constant (posteriors, energy-based models), it is consistent with an $O(1/n)$ bias, and the same weights estimate $Z$.
4. For $p = \mathcal N(0, 1)$ and $q = \mathcal N(0, s^2)$, for which $s$ is the IS variance of the estimate of $\mathbb{E}[X^2]$ finite?
   Answer

   The second moment is $\int x^4 p^2/q\,dx$ with $p^2/q \propto \exp\big(-x^2(1 - \frac{1}{2s^2})\big)$. This is integrable (also times $x^4$) exactly when $1 - \frac{1}{2s^2} > 0$, that is $s^2 > 1/2$.
5. For which of $h(u) = u$, $h(u) = e^u$, $h(u) = (u - \frac12)^2$ do antithetic pairs help?
   Answer

   $h(u) = u$: $h(u) + h(1 - u) = 1$, so the variance is 0. $e^u$: monotone, so $\rho < 0$ and it helps a lot ($\rho = -0.968$). $(u - \frac12)^2$: symmetric, $\rho = 1$, so the variance per evaluation doubles.
6. A control variate has correlation $\rho = 0.9$ with $f$. How many times fewer draws do you need? And if $\rho = -0.9$?
   Answer

   The variance factor is $1 - \rho^2 = 0.19$, so about $1/0.19 \approx 5.3$ times fewer draws. The sign does not matter: $\beta^*$ takes the sign of the covariance.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the standard error of the plain Monte Carlo estimate of $\mathbb{E}[f(X)]$ from $n$ i.i.d. draws?

   Show answer

2. What is the relative standard error of plain Monte Carlo for an event of probability $\mu$?

   Show answer

3. In rejection sampling, when do you accept a proposal $x \sim q$, and what is the acceptance rate?

   Show answer

4. Write the importance sampling estimator of $\mathbb{E}_p[f(X)]$.

   Show answer

5. Why must an importance sampling proposal have heavier tails than the target?

   Show answer

6. What is the self-normalised IS estimator, and is it unbiased?

   Show answer

7. Define the effective sample size of importance weights.

   Show answer

8. What is the optimal coefficient of a control variate, and the resulting variance factor?

   Show answer

9. When do antithetic variables reduce the variance, and when do they increase it?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Estimate $\pi$ with a random number generator. How accurate is it?" — Fraction of points in the disk times 4; standard error $4\sqrt{(\pi/4)(1 - \pi/4)/n} \approx 1.64/\sqrt n$.
• "You have a density up to a constant. How do you sample it?" — Inverse CDF in 1-D if you can integrate and invert; rejection with an envelope in low dimension; otherwise MCMC (lesson 15), or SNIS if you only need expectations.
• "Why does importance sampling fail in high dimension?" — The variance of the weights grows exponentially with the dimension (ESS $= n\,e^{-\|m\|^2}$ for a mean shift $m$), so one weight dominates.
• "Why is speculative decoding exact?" — Accept with $\min(1, p/q)$ and, on rejection, sample from the normalised residual $\max(0, p - q)$; the two parts add up to $p$.
• "Importance sampling or MCMC for a posterior?" — IS from a Laplace or variational approximation with heavier tails when the dimension is low and a good proposal exists: independent draws, the evidence $Z$ and a weight diagnostic. MCMC when the posterior is high-dimensional or has a shape no simple proposal covers.
• "What is a baseline in policy gradients?" — A control variate: subtracting $b$ keeps the estimator unbiased because $\mathbb{E}[\nabla\log\pi] = 0$, and reduces the variance when $b$ is close to the optimal value; a bad $b$ can increase it.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. log_sum_exp and effective_sample_size from log-weights, stable for weights like $e^{\pm 1000}$.
2. importance_sampling for $P(X > 4)$ with its standard error, compared with scipy.stats.norm.sf.
3. self_normalized_is on the unnormalised mixture $3\,\mathcal N(-2, 1) + \mathcal N(3, 0.5^2)$ (so $Z = 4$): the mean, and $\hat Z$.
4. A rejection sampler for the half-normal with an Exp(1) proposal: the acceptance rate against the optimal $M$.
5. control_variate_mean and antithetic_mean, with their variance-reduction factors for $\ln 2$.

11. Watch and read #

Importance SamplingMutual Information · 13 min

Why sampling from a different distribution can reduce variance, and how the weights correct the bias.

Monte Carlo Simulations : Data Science Basicsritvikmath · 19 min

The basic idea of estimating quantities by random sampling, with simple examples.

Machine learning - Importance sampling and MCMC INando de Freitas · 1 h 16 min

Monte Carlo integration, importance sampling and the move to MCMC, from a machine-learning course.

• Art Owen, Monte Carlo theory, methods and examples (free chapters): ch. 2 (simple Monte Carlo), ch. 4 §4.1–4.4 (inversion) and §4.7 (acceptance–rejection), ch. 8 §8.2 (antithetics) and §8.9 (control variates), ch. 9 §9.1–9.3 (importance sampling, self-normalisation, the ESS).
• Kevin Murphy, Probabilistic Machine Learning: Advanced Topics (free PDF): ch. 11 (Monte Carlo methods: inverse CDF, rejection sampling, importance sampling, variance reduction).
• Yaniv Leviathan, Matan Kalman and Yossi Matias, Fast Inference from Transformers via Speculative Decoding (2023): rejection sampling at the scale of LLMs.
• Yuri Burda, Roger Grosse and Ruslan Salakhutdinov, Importance Weighted Autoencoders (2016): importance sampling inside a variational bound.