Sampling and Monte Carlo Methods

Many quantities in machine learning are expectations: the expected loss, the posterior predictive, the normalising constant of a distribution, the value of a policy. Most of these integrals have no closed form. Monte Carlo methods estimate them by drawing random samples. The idea is old — it was named during the Manhattan Project after the casino in Monaco — and it is one of the most useful ideas in computational science.

Monte Carlo estimation#

To estimate $\mu = \mathbb{E}_{p}[f(X)] = \int f(x)p(x)\,dx$, draw $x_1, \dots, x_N \sim p$ and average:

\hat\mu_N = \frac{1}{N}\sum_{i=1}^{N} f(x_i)

The estimator is unbiased, and its standard error is $\sigma_f/\sqrt{N}$. Crucially, this rate does not depend on dimension — unlike grid-based numerical integration, whose cost explodes exponentially with dimension. That is why Monte Carlo dominates in high dimensions.

python

import numpy as np

rng = np.random.default_rng(0)
# Estimate pi: fraction of random points in the unit square inside the quarter circle
for N in [100, 10_000, 1_000_000]:
    pts = rng.random((N, 2))
    est = 4 * np.mean((pts**2).sum(1) <= 1)
    print(f"N={N:>9}: pi ~ {est:.5f}   (error {abs(est - np.pi):.5f})")

Every factor of 100 in samples gives one more correct digit — the $1/\sqrt{N}$ law in action.

Generating samples#

Inverse transform sampling#

If $F$ is the CDF of $X$ and $U \sim \text{Uniform}(0,1)$, then $F^{-1}(U)$ has distribution $F$. For the exponential distribution, $F^{-1}(u) = -\ln(1 - u)/\lambda$.

Rejection sampling#

To sample from a target $p(x)$ known up to a constant, use a proposal $q(x)$ with $M q(x) \ge \tilde{p}(x)$ everywhere. Draw $x \sim q$ and $u \sim U(0,1)$; accept if $u < \tilde{p}(x)/(Mq(x))$. It is exact but becomes hopelessly inefficient in high dimensions because the acceptance rate collapses.

Importance sampling#

Sometimes we cannot sample from $p$, or sampling from $p$ rarely hits the region that matters (rare events). Sample from a proposal $q$ instead and reweight:

\mathbb{E}_p[f(X)] = \mathbb{E}_q\left[f(X)\frac{p(X)}{q(X)}\right] \approx \frac{1}{N}\sum_{i=1}^{N} f(x_i)\,w_i, \qquad w_i = \frac{p(x_i)}{q(x_i)}

When $p$ is only known up to a constant, normalise the weights (self-normalised importance sampling). The variance depends heavily on how well $q$ matches $|f|p$; a poor proposal produces a few enormous weights and a useless estimate. The effective sample size $\text{ESS} = (\sum w_i)^2/\sum w_i^2$ diagnoses this.

Markov chain Monte Carlo (MCMC)#

For complex high-dimensional distributions — such as Bayesian posteriors known only up to a normalising constant — we construct a Markov chain whose stationary distribution is the target. After a burn-in period, the chain's states are (correlated) samples from $p$.

Metropolis–Hastings#

From the current state $x$, propose $x' \sim q(x' \mid x)$ and accept with probability

A = \min\left(1,\; \frac{\tilde{p}(x')\,q(x \mid x')}{\tilde{p}(x)\,q(x' \mid x)}\right)

The normalising constant cancels — we only need the unnormalised density $\tilde{p}$. This is why MCMC is so powerful for Bayesian inference, where the evidence $p(\mathcal{D})$ is intractable.

python

import numpy as np

def log_target(x):          # unnormalised mixture of two Gaussians
    return np.logaddexp(-0.5 * (x + 2) ** 2, -0.5 * ((x - 3) / 0.7) ** 2)

rng = np.random.default_rng(1)
x, samples, accepted = 0.0, [], 0
for t in range(60_000):
    prop = x + rng.normal(0, 1.5)                  # symmetric proposal: q terms cancel
    if np.log(rng.random()) < log_target(prop) - log_target(x):
        x, accepted = prop, accepted + 1
    if t >= 10_000:                                # discard burn-in
        samples.append(x)
print("acceptance rate:", accepted / 60_000, " mean:", np.mean(samples).round(3))

Gibbs sampling#

Update one variable at a time from its full conditional $p(x_j \mid \mathbf{x}_{-j})$. It is a special case of Metropolis–Hastings with acceptance probability one, ideal when conditionals are easy (e.g. conjugate models, LDA topic models).

Hamiltonian Monte Carlo#

Random-walk proposals explore high-dimensional spaces slowly. HMC uses gradients of $\log p$ to simulate physical dynamics and make long, informed moves. The NUTS variant tunes itself automatically and powers probabilistic programming tools such as Stan, PyMC and NumPyro.

Diagnosing MCMC#

Trace plots should look like "fuzzy caterpillars", not slow drifts.
Run multiple chains from different starts and compare with $\hat{R}$ (should be close to 1.0).
Report effective sample size, since samples are autocorrelated.

Monte Carlo in deep learning#

SGD itself is a Monte Carlo estimate of the full gradient.
Dropout at test time (MC dropout) approximates Bayesian model averaging.
VAEs estimate the ELBO with sampled latent variables.
Policy gradients estimate expected returns from sampled trajectories.
Diffusion models generate data by iteratively sampling from learned conditional distributions.

#sampling#monte carlo#importance sampling#mcmc#metropolis-hastings

Written by

Janin A Apurba

B.Sc. in CSE, AUST · Advanced ICT Officer, CNRS-UNHCR. Teaching AI, ML and Deep Learning to the next generation of engineers and researchers.