Sampling methods

The last lesson ended on a frustration: NumPy samples a Normal or a Binomial in one line, but the distribution you most want in Bayesian work — the posterior you just learned to build — has no name and no built-in sampler. This lesson is the way around that.

Two questions come up over and over in ML, and both are answered by sampling. One: given a distribution, draw random examples from it — even an unnamed one. Two: estimate an expectation E[f(X)] — an average, a probability, an integral — that you cannot compute in closed form. When the integral is intractable (which, in high dimensions, is almost always), you stop trying to solve it and start sampling it.

Monte Carlo: estimate by averaging samples

The core idea is Monte Carlo estimation: to approximate E[f(X)], draw N samples x_1 … x_N and average:

E[f(X)]  ≈  (1/N) · ( f(x_1) + f(x_2) + ... + f(x_N) )

The error shrinks like 1/√N — and, remarkably, that rate does not depend on the dimension. A grid-based integral needs exponentially more points as dimensions grow (the curse of dimensionality). Monte Carlo’s 1/√N is the same in 2 dimensions or 2,000, which is exactly why it dominates high-dimensional ML. The classic demonstration is estimating π by throwing darts at a square and counting how many land in the quarter circle:

import numpy as np

rng = np.random.default_rng(0)
for N in [100, 1000, 10000, 100000]:
    pts = rng.random((N, 2))
    inside = int(((pts**2).sum(axis=1) <= 1).sum())
    print(f"N={N:>6}: pi ~ {4 * inside / N:.4f}")

N=   100: pi ~ 2.9200
N=  1000: pi ~ 3.2920
N= 10000: pi ~ 3.1460
N=100000: pi ~ 3.1431

At N = 100 the estimate is a noisy 2.92; by N = 100,000 it’s 3.143. Tenfold more samples roughly halves the error (1/√N) — slow but utterly general.

Getting the samples

Monte Carlo assumes you can draw from the distribution. When there’s no built-in sampler, a few techniques build one:

• Inverse-transform. If you can invert the CDF F, sample u ~ Uniform(0,1) and return F⁻¹(u). Exact and cheap — but only when the inverse CDF is available (mostly 1-D, known distributions).
• Rejection sampling. Draw from an easy proposal that envelopes the target, then accept each draw with probability target/envelope (reject otherwise). Simple, but it wastes most samples when the envelope fits poorly — and the waste explodes in high dimensions.
• Importance sampling. You don’t even have to sample the target. Sample from a convenient proposal q, then reweight each sample by w = p(x)/q(x) to estimate expectations under p. It’s how off-policy reinforcement learning reuses data from a different policy.

MCMC: when the distribution is intractable

The hard, common case in ML is a distribution you can only evaluate up to a constant — a Bayesian posterior p(θ | data) ∝ p(data | θ) p(θ), whose normalizer is an intractable integral.

Markov Chain Monte Carlo (MCMC) handles it: build a Markov chain whose stationary distribution is the target, then walk the chain and collect its states as samples. Metropolis-Hastings is the canonical recipe — propose a move, accept it with a probability that depends only on the ratio of densities, so the unknown normalizer cancels. Over time, the walk visits states in proportion to the target. MCMC is the workhorse of Bayesian inference precisely because it never needs the normalizing constant.

In one breath

• Two needs: draw samples from a distribution, or estimate an expectation E[f(X)] you can’t integrate.
• Monte Carlo: E[f(X)] ≈ (1/N) Σ f(x_i), error ~1/√N — dimension-independent, which is why it beats grids in high-D (the π-by-darts demo: 2.92 → 3.143 as N grows).
• To get samples: inverse-transform (invert the CDF, exact, mostly 1-D), rejection (proposal + accept/reject, wasteful in high-D), importance sampling (sample q, reweight by p/q).
• MCMC (Metropolis-Hastings) samples intractable, unnormalized distributions (Bayesian posteriors) by walking a Markov chain whose stationary distribution is the target — the normalizer cancels in the acceptance ratio.
• It powers Bayesian inference, diffusion, RL, dropout, VAEs, and the bootstrap.

Practice

A question to carry forward

Look back at the trick that made MCMC work — the one we waved at without unpacking. To sample an intractable posterior, we built a Markov chain whose stationary distribution is the target, then walked it. Three loaded phrases, all taken on faith: what is a Markov chain, what does it mean for a distribution to be “stationary,” and why on earth should aimlessly walking a chain leave you holding samples from a specific distribution?

So here is the thread onward. A random variable is one uncertain number; but MCMC needed randomness that unfolds over time — a whole sequence of linked random states. What is a stochastic process, and what is the Markov property that makes the next step depend only on the current one? Why does such a chain settle into a single stationary distribution no matter where it began — the very fact that quietly turns “take a random walk” into “draw a sample”?