Monte Carlo methods

The law of large numbers turned every expectation into a limit of averages: draw samples, average, and the error falls like $1/\sqrt n$. Two questions remain before that recipe is a method. How are samples from a given distribution produced at all, when a computer only supplies uniform random numbers? And what can be done when $1/\sqrt n$ is too slow, as it is for rare events? This note answers both: inversion, transformation and rejection produce samples; importance sampling and variance reduction make each sample count for more.

Notation. The target is $\mu = E_p[g(X)]$ for a density $p$; $\hat\mu_n = \frac1n \sum g(X_i)$ with $X_i \sim p$. A proposal $q$ is a density that is easy to sample.

Plain Monte Carlo and its error#

With $X_i$ independent from $p$ and $E_p\lvert g(X) \rvert < \infty$, $\hat\mu_n$ is unbiased and converges to $\mu$. If moreover $\sigma^2 = \operatorname{Var}_p(g(X))$ is finite, by the central limit theorem it is approximately $N(\mu, \sigma^2/n)$. A Monte Carlo estimate should always be reported with its error bar, $\hat\mu_n \pm 1.96\, \hat\sigma/\sqrt n$, computed from the same samples. The error depends on $\sigma$ and $n$ and not directly on the dimension of $X$, which is why Monte Carlo is the method of choice for high-dimensional integrals. Without a finite variance the error bar is meaningless, even when the estimate converges: heavy-tailed $g(X)$, and importance weights with infinite variance below, need other diagnostics or robust estimators such as the median of means. Everything below either produces the samples or reduces $\sigma$.

Producing samples#

A random number generator supplies independent $U(0, 1)$ numbers. Three techniques turn them into samples from other distributions.

• Inversion. $F^{-1}(U)$ has CDF $F$, as shown in the functions of random variables note: $-\ln(1 - U)$ is exponential, and a lookup in cumulative sums samples any discrete distribution.
• Transformation. When $F^{-1}$ has no formula, another map may. The Box–Muller transform takes two independent uniforms to a standard normal: $Z = \sqrt{-2\ln(1 - U_1)} \cos(2\pi U_2)$. It works because the pair $(Z_1, Z_2)$ of independent normals has a rotation-invariant density, so its angle is uniform and independent of the radius, and its squared radius, a sum of two squared standard normals, is exponential with mean 2.
• Rejection. For a density known only up to a constant, or with no convenient sampling formula, propose and accept.

Rejection sampling#

Choose a proposal $q$ that is easy to sample and a constant $M$ with $p(x) \le M q(x)$ everywhere. Draw $X \sim q$ and $U \sim U(0, 1)$, and accept $X$ if $U M q(X) \le p(X)$; otherwise discard it and try again. The point $(X, U M q(X))$ is uniform under the envelope $M q$, so the accepted points are uniform under $p$, and their $x$-coordinates have density $p$:

Figure

Rejection sampling

The method needs $p$ only up to a constant factor, since a constant can be absorbed into $M$. Its cost is set by the acceptance rate, $1/M$ for a normalized target: on average $M$ proposals per accepted sample. That is $2.46$ for the example, but it typically grows exponentially with the dimension, or faster. Sampling the unit ball by rejection from the surrounding cube accepts 79% of proposals in two dimensions, 0.25% in ten and $2.5 \cdot 10^{-8}$ in twenty: almost all the volume of a high-dimensional cube lies in its corners, as the multivariate Gaussian note found for the shell of a Gaussian. Rejection sampling is a tool for low dimensions and for the inner loops of other samplers.

Importance sampling#

Plain Monte Carlo spends its samples where $p$ is large, which may be far from where $g$ matters. Importance sampling draws from a proposal $q$ instead and corrects with weights:

$$E_p[g(X)] = \int g(x) \frac{p(x)}{q(x)}\, q(x)\, dx = E_q\big[g(X)\, w(X)\big], \qquad w(x) = \frac{p(x)}{q(x)} ,$$

valid whenever $q > 0$ wherever $g p \ne 0$. The estimator $\frac1n \sum g(X_i) w(X_i)$ with $X_i \sim q$ is unbiased, and its variance depends on how well $q$ matches $\lvert g \rvert p$; the ideal proposal is proportional to it.

The gain is largest for rare events. To estimate $P(Z > 4) = 3.17 \cdot 10^{-5}$ for a standard normal $Z$, a million plain samples see only 28 events, against the $31.7$ expected: a relative standard error of about 18%, and a realized error of 12% in this run. Sampling from $q = N(4, 1)$ puts half the draws in the event, each weighted by $p/q = e^{8 - 4x}$:

Figure

Importance sampling for a rare event

Ten thousand draws from $q$ estimate $3.27 \cdot 10^{-5}$ with a relative standard error of 2%, a hundred times fewer samples for eight times less error; the realized error of this estimate, known here because the exact value is, is 3.3%. The same weights are useless for other questions: their effective sample size $(\sum w_i)^2 / \sum w_i^2$ is $5.1$ out of 10,000, so estimating, say, the mean of $Z$ with them would be hopeless. A proposal is good for a specific integrand, and a proposal with lighter tails than $p$ can have infinite variance.

When $p$ is known only up to a constant, as a posterior is, the weights are normalized by their sum (self-normalized importance sampling); the estimate is then slightly biased but consistent. Importance weights appear throughout machine learning: off-policy evaluation in reinforcement learning reweights logged actions by the ratio of the new policy to the old, and the importance-weighted bound of variational autoencoders averages the importance weights of several samples inside the logarithm.

Variance reduction#

Even with the right distribution, the variance of the estimator can often be cut by using what is already known.

• Control variates. If a quantity $h(X)$ with known mean is correlated with $g(X)$, estimate $E[g(X) + c\,(h(X) - E h(X))]$ instead, with $c = -\operatorname{Cov}(g, h)/\operatorname{Var}(h)$. To compute $E[e^U] = e - 1$ for $U \sim U(0, 1)$, the control $h(U) = U$ with mean $1/2$ cuts the variance per sample from $0.2446$ to $0.0040$, a factor of 60, for free. In policy-gradient reinforcement learning, subtracting a baseline from the reward is a control variate for the gradient estimate, and it is one of the main reasons such methods work in practice.
• Antithetic pairs. Use $U$ and $1 - U$ together: for monotone $g$ the two values are negatively correlated, and their average has lower variance than two independent samples.
• Stratification and quasi-Monte Carlo. Spread the samples evenly over the space instead of at random; in low dimensions, low-discrepancy sequences reach errors close to $1/n$ instead of $1/\sqrt n$.

Gradient estimation in machine learning is the same problem. The gradient of an expectation $E_{q_\theta}[f(Z)]$ can be estimated by the score-function (REINFORCE) estimator, unbiased but high-variance, or by the reparametrization trick of the variational autoencoder, which writes $Z = \mu_\theta + \sigma_\theta \varepsilon$ and usually has far lower variance.

What changed#

	Plain Monte Carlo	Monte Carlo methods
Samples from $p$	assumed available	inversion, transformation, rejection
Error	$\sigma/\sqrt n$, reported with the estimate	the same rate, smaller $\sigma$
Rare events	most samples wasted	importance sampling from a proposal near the event
Known structure	ignored	control variates, antithetic pairs, stratification
High dimensions	the rate does not depend on $d$	but rejection rates and importance weights degrade with $d$

Limitations#

• Good proposals are hard to find in high dimensions. Both rejection sampling and importance sampling need $q$ close to $p$; in many dimensions a small mismatch in every coordinate compounds, and acceptance rates and effective sample sizes collapse exponentially.
• Importance sampling can fail silently. A proposal with lighter tails than the target can give estimates whose variance is infinite: they look stable for a long time and then jump.
• Exact samples from complicated distributions are rarely available. Posteriors of real models, energy-based models and distributions over structures need a different approach, the subject of the next note.

What you should remember#

• Report a Monte Carlo estimate with its error bar $\pm 1.96\,\hat\sigma/\sqrt n$, valid for independent draws when $g(X)$ has a finite variance; the rate does not depend on the dimension.
• Samples come from inversion ($F^{-1}(U)$), transformations (Box–Muller) or rejection (accept $X \sim q$ with probability $p(X)/(Mq(X))$, acceptance rate $1/M$).
• Rejection sampling needs $p$ only up to a constant but its acceptance rate collapses in high dimensions.
• Importance sampling: $E_p[g] = E_q[g\, p/q]$; it can make rare events cheap; check the effective sample size $(\sum w)^2/\sum w^2$.
• Control variates subtract a correlated quantity with known mean; baselines in policy gradients are control variates.

Code#

Rejection sampling, Box–Muller, a rare event by importance sampling and a control variate:

import numpy as np, math
from scipy import stats

rng = np.random.default_rng(0)

# rejection sampling from Beta(2, 5) with a uniform proposal and M = max p = p(0.2)
p = stats.beta(2, 5).pdf
M = p(0.2)
x, u = rng.random(100_000), rng.random(100_000)
sample = x[u * M <= p(x)]
print(round(float(M), 4), round(len(sample) / 100_000, 3), round(float(sample.mean()), 3), round(2 / 7, 3))   # 2.4576 0.408 0.286 0.286

# the share of a cube taken by the inscribed ball: the acceptance rate of rejection sampling in d dimensions
print([f'{math.pi ** (d / 2) / math.gamma(d / 2 + 1) / 2 ** d:.2g}' for d in (2, 5, 10, 20)])   # ['0.79', '0.16', '0.0025', '2.5e-08']

# Box-Muller: two uniforms give a standard normal
u1, u2 = rng.random((2, 500_000))
z = np.sqrt(-2 * np.log(1 - u1)) * np.cos(2 * np.pi * u2)
print(round(float(z.mean()), 3), round(float(z.std()), 3))           # 0.0 1.003

# a rare event, P(Z > 4): a million plain samples against 10,000 importance samples from N(4, 1)
true = stats.norm.sf(4)
plain = rng.normal(size=1_000_000) > 4
y = rng.normal(4, 1, 10_000)
w = np.exp(8 - 4 * y)                                                # p(y) / q(y)
terms = (y > 4) * w
print(f'{true:.3g}', int(plain.sum()), f'{terms.mean():.3g}', round(float(terms.std() / math.sqrt(1e4) / true), 3))   # 3.17e-05 28 3.27e-05 0.022
print(round(float(w.sum() ** 2 / (w ** 2).sum()), 1))               # 5.1: the effective sample size of these weights

# a control variate for E[e^U] = e - 1: the variance per sample without and with it
U = rng.random(10_000)
Y = np.exp(U)
c = -np.cov(Y, U)[0, 1] / U.var(ddof=1)
print(round(float(Y.mean()), 4), round(float(Y.var()), 4), round(float((Y + c * (U - 0.5)).var()), 4), round(math.e - 1, 4))   # 1.7181 0.2446 0.004 1.7183

Next#

Rejection and importance sampling need a proposal close to the target, and in high dimensions there usually is none. The last note of the course builds instead a Markov chain whose stationary distribution is the target, using only the unnormalized density: the Metropolis–Hastings algorithm, Gibbs sampling, and Langevin dynamics, the bridge to diffusion models.

References#

• C. P. Robert, G. Casella. Monte Carlo Statistical Methods, 2nd ed. Springer, 2004. Chapter 2 (random variable generation, rejection), Chapter 3 (Monte Carlo integration, importance sampling), Chapter 4 (variance reduction).
• A. B. Owen. Monte Carlo Theory, Methods and Examples. 2013. Chapters 4, 8 and 9 (non-uniform random variables, variance reduction, importance sampling).
• G. E. P. Box, M. E. Muller. A Note on the Generation of Random Normal Deviates. The Annals of Mathematical Statistics, 29(2), 1958.
• J. von Neumann. Various Techniques Used in Connection with Random Digits. National Bureau of Standards Applied Mathematics Series, 12, 1951. The rejection method.
• R. J. Williams. Simple Statistical Gradient-Following Algorithms for Connectionist Reinforcement Learning. Machine Learning, 8, 1992. REINFORCE and baselines.