Sampling and Monte Carlo methods

Drawing from distributions by inverse CDF and rejection, estimating expectations with importance sampling, and reducing variance with control variates.

Sampling from a distribution

If you can invert the cumulative distribution function, you can sample: draw u uniform on [0, 1] and return F^-1(u). This is how exponential, categorical and many discrete distributions are sampled in practice.

import numpy as np
rng = np.random.default_rng(0)

# inverse CDF: exponential with rate lambda
def sample_exponential(lam, size):
    u = rng.random(size)
    return -np.log(1 - u) / lam

x = sample_exponential(2.0, 100_000)
x.mean(), 1 / 2.0                        # 0.5003, 0.5  (mean is 1/lambda)

# Box-Muller: two uniforms become two independent Gaussians
def sample_normal(size):
    u1 = rng.random(size)
    u2 = rng.random(size)
    return np.sqrt(-2 * np.log(u1)) * np.cos(2 * np.pi * u2)

sample_normal(100_000).std()             # about 1.0

# categorical: one pass over the cumulative sum
weights = np.array([0.5, 0.3, 0.2])
draws = rng.choice(len(weights), size=10_000, p=weights)
np.bincount(draws) / 10_000

• Antithetic pairs, u and 1-u, cancel much of the sampling noise for monotone functions at no extra cost in distributional correctness.
• The reparameterisation trick (mu + sigma * eps with eps standard normal) makes a sample differentiable with respect to the distribution's parameters — this is what makes VAEs trainable.
• Sampling is O(1) per draw in expectation for the methods above; a naive rejection scheme can be arbitrarily slow when the proposal fits badly.
• Reproducibility requires a seeded generator object (np.random.default_rng(seed)), not the global legacy state.

Rejection and importance sampling

# rejection sampling from a target we know only up to a constant
def target(x):                       # unnormalised: a mixture of two bumps
    return np.exp(-0.5 * ((x - 2) / 0.7) ** 2) + 0.5 * np.exp(-0.5 * ((x + 2) / 1.0) ** 2)

def rejection_sample(size, proposal_scale=4.0, M=1.2):
    out = []
    while len(out) < size:
        x = rng.uniform(-proposal_scale, proposal_scale)   # proposal / (2*scale)
        if rng.uniform(0, M * 1.0 / (2 * proposal_scale)) < target(x):
            out.append(x)
    return np.array(out)

samples = rejection_sample(20_000)
samples.mean()                                       # about 0.4 (mean of each bump)

Importance sampling does not sample from the target at all. It draws from a convenient proposal q and reweights every sample by p(x)/q(x), giving an unbiased estimate of an expectation under p. It is the standard tool for estimating a rare-event probability whose distribution you cannot sample directly.

# estimate E_p[f] using samples from a shifted proposal
f = lambda x: (x > 3).astype(float)         # a rare event under p

xs_chosen = rng.normal(0, 1, 200_000)
xs_prop = rng.normal(3, 1, 200_000)         # proposal centred on the event

def p_density(x):  return np.exp(-0.5 * x ** 2) / np.sqrt(2 * np.pi)
def q_density(x):  return np.exp(-0.5 * (x - 3) ** 2) / np.sqrt(2 * np.pi)

def estimate(samples, qd):
    w = p_density(samples) / qd(samples)
    return (w * f(samples)).mean(), w.std() / np.sqrt(len(w))

print("direct ", estimate(xs_chosen, p_density))
print("proposal", estimate(xs_prop, q_density))

Monte Carlo error and variance reduction

# a control variate: subtract a function with a known mean
def control_variate(y, x, known_mean):
    cov = np.cov(y, x)[0, 1]
    c = cov / np.var(x)
    return y - c * (x - known_mean)

# estimator of pi by counting points in a quarter circle
u = rng.random((1_000_000, 2))
inside = ((u ** 2).sum(axis=1) <= 1.0).astype(float)
pi_hat = 4 * inside.mean()
stderr = 4 * inside.std() / np.sqrt(len(inside))
pi_hat, stderr, np.pi

• Monte Carlo error falls as sigma / sqrt(n) regardless of dimension. That dimension-independence is why the method is used for high-dimensional integrals that no quadrature rule can handle.
• To halve the standard error you need four times the samples: 100x more compute buys 10x more precision, a poor exchange once each sample is expensive.
• Antithetic variates and control variates reduce sigma itself and are effectively free when the structure fits.
• For a model's predictive uncertainty, the practical trick is to run the forward pass several times with dropout active (Monte Carlo dropout) and look at the spread of the outputs.

FAQ

Related

Probability and distributions Statistical estimation

Last refreshed 2026-09-18.