Monte Carlo and sampling basics

"When I cannot solve the integral, I ask ten thousand dice to vote on the answer, and I trust the average of their opinions."

A Statistically Enfranchised AI Agent

This section assumes you are comfortable with random variables, distributions, and expectations from Chapter 4, sections 4.1 through 4.3, and it uses the aleatoric-versus-epistemic distinction of section 4.4 to decide what to sample. If you have not met the noise models that produce these distributions, section 2 of Chapter 2 supplies them. No measure theory is required; a working knowledge of the law of large numbers is enough.

The Monte Carlo estimator: averaging your way to an integral

What. Monte Carlo estimation answers the question "what is \(\mathbb{E}_{x\sim p}[f(x)]\)?" by drawing \(N\) independent samples \(x_1,\dots,x_N\) from \(p\) and computing the sample mean

Why it works. The estimator \(\hat{\mu}_N\) is unbiased: its expectation equals the true integral for every \(N\). By the law of large numbers it converges to the true value as \(N\) grows, and by the central limit theorem its error is approximately Gaussian with standard deviation \(\sigma/\sqrt{N}\), where \(\sigma^2 = \operatorname{Var}_p[f(x)]\). That \(\sqrt{N}\) is the single most important fact in this section. To halve your error you must quadruple your samples. The error does not care how many dimensions \(x\) has, which is why Monte Carlo beats grid-based numerical integration the moment your state vector grows past three or four dimensions. A 12-dimensional IMU-plus-bias state would need a grid of astronomical size; the same accuracy from Monte Carlo needs only enough samples to tame the variance.

How you use it. Nearly every quantity you want is an expectation in disguise. A probability \(P(A)\) is \(\mathbb{E}[\mathbb{1}_A]\), the expected value of an indicator. A quantile is found by sampling and sorting. A predictive mean is \(\mathbb{E}[f(x)]\); a predictive variance is \(\mathbb{E}[f(x)^2]-\mathbb{E}[f(x)]^2\). Choose \(f\), draw from \(p\), average.

Turning uniform noise into any distribution

What. Your hardware random-number generator gives you one primitive: uniform draws on \([0,1)\). Everything else is built from that primitive by transformation.

How, method one: inverse transform. If \(F\) is the cumulative distribution function of a scalar random variable and \(U\sim\text{Uniform}(0,1)\), then \(X = F^{-1}(U)\) has exactly the distribution you wanted. This is why an exponential inter-arrival time for a Poisson event stream is generated as \(-\lambda^{-1}\ln(1-U)\): the log is the inverse CDF. Inverse transform is exact and cheap whenever you can invert \(F\).

How, method two: rejection sampling. When you cannot invert \(F\) but you can evaluate an unnormalized density \(\tilde{p}(x)\), wrap it in a proposal \(q(x)\) you can sample, scaled so that \(M\,q(x)\ge\tilde{p}(x)\) everywhere. Draw \(x\sim q\), accept it with probability \(\tilde{p}(x)/(M\,q(x))\), and repeat. Accepted samples are distributed exactly as \(p\). The catch is efficiency: the acceptance rate is \(1/M\), so a loose envelope in high dimensions wastes almost every draw.

How, method three: importance sampling. Sometimes you cannot or should not draw from \(p\) at all. Instead draw from a convenient \(q\) and reweight: \(\mathbb{E}_p[f] = \mathbb{E}_q[f(x)\,w(x)]\) with weights \(w(x)=p(x)/q(x)\). This is the mathematical heart of the particle filter you will build in Chapter 10, where \(q\) is the motion model and the weights carry the measurement likelihood. Importance sampling shines when you deliberately over-sample a rare but critical region, such as the tail where a fall-detection alarm should fire, then correct the bias with weights.

How many samples? Error bars on your error bars

What. A Monte Carlo estimate without an uncertainty is half an answer. Because \(\hat{\mu}_N\) is itself a random variable, you can estimate its standard error from the very same samples: \(\widehat{\text{SE}} = \hat{\sigma}/\sqrt{N}\), where \(\hat{\sigma}\) is the sample standard deviation of the \(f(x_i)\) values. Report \(\hat{\mu}_N \pm 1.96\,\widehat{\text{SE}}\) for an approximate 95% interval.

Why it matters for sensing. When you estimate a rare-event probability such as a \(10^{-4}\) false-alarm rate, naive sampling needs on the order of \(10^6\) draws just to see a handful of events, and the relative error stays enormous until you do. This is exactly where importance sampling earns its keep: bias the proposal toward the tail, and the same accuracy arrives with orders of magnitude fewer draws. The estimator that turns raw draws into an answer with a defensible interval is the code below.

import numpy as np

rng = np.random.default_rng(42)

def mc_estimate(f, sampler, n, batch=100_000):
    """Streaming Monte Carlo mean with a 95% standard-error interval.
    Batches keep memory flat, which matters on edge hardware."""
    total, total_sq, count = 0.0, 0.0, 0
    while count < n:
        m = min(batch, n - count)
        vals = f(sampler(m))          # f applied to a batch of draws
        total    += vals.sum()
        total_sq += (vals ** 2).sum()
        count    += m
    mean = total / count
    var  = total_sq / count - mean ** 2
    se   = np.sqrt(max(var, 0.0) / count)
    return mean, 1.96 * se

# Estimate P(||accel bias|| > 0.05 g) for a 3-axis IMU with correlated bias.
cov = np.array([[4e-4, 1e-4, 0.0],
                [1e-4, 4e-4, 0.0],
                [0.0,  0.0,  9e-4]])   # units of g^2
L = np.linalg.cholesky(cov)
sample_bias = lambda m: (L @ rng.standard_normal((3, m))).T
exceeds = lambda b: (np.linalg.norm(b, axis=1) > 0.05).astype(float)

p_hat, halfwidth = mc_estimate(exceeds, sample_bias, n=2_000_000)
print(f"P(exceed) = {p_hat:.4f} +/- {halfwidth:.4f}")

The example applies the estimator to a real specification question: how often does an IMU's correlated bias vector exceed a magnitude limit? The Cholesky factor turns independent Gaussian draws into correlated ones, a trick you will reuse constantly when the noise channels of a sensor are coupled. The returned half-width tells you whether two million samples were enough or whether you should keep drawing.

Variance reduction and low-discrepancy sampling

What. Since error scales as \(\sigma/\sqrt{N}\) and buying more \(N\) is expensive on an edge device, the smart move is to shrink \(\sigma\). Three techniques recur throughout sensing.

Antithetic variates. For every draw \(x\), also use its mirror image (for a symmetric distribution, \(-x\) after centering). Positively correlated errors cancel, and a symmetric integrand can see its variance drop sharply at zero extra sampling cost. Control variates subtract a correlated quantity whose expectation you know analytically, replacing part of the noisy estimate with an exact term. Stratified sampling partitions the domain and samples each stratum in proportion, guaranteeing coverage rather than trusting luck.

Quasi-Monte Carlo. Instead of pseudo-random points, use a deterministic low-discrepancy sequence such as Sobol or Halton that fills space more evenly than random draws ever do. For smooth integrands the error can improve toward \(1/N\) rather than \(1/\sqrt{N}\), a decisive speedup. This is the workhorse behind the digital-twin and synthetic-data sweeps of Chapter 55, where you want to cover a parameter space of sensor placements or lighting conditions without wasting simulation budget on clumped points.

When to reach for each. Use antithetic and control variates when you have structure you can exploit and evaluations are cheap. Use quasi-Monte Carlo for smooth, moderate-dimensional integrals such as calibration sweeps. Fall back to plain pseudo-random sampling when the integrand is rough, discontinuous, or when you need the honest, well-understood \(\sqrt{N}\) error bars that the central limit theorem guarantees; low-discrepancy sequences complicate confidence-interval construction because their points are not independent.

What's Next

In Chapter 5, we leave the world of clean synthetic draws and confront real recorded data: how to build datasets from sensor streams without leaking future information into the past, how to split by device and by time rather than at random, and how the sampling intuition from this section informs honest train-test protocols. The same law of large numbers that made our estimators trustworthy will now warn us about what happens when samples are correlated and not the independent draws Monte Carlo assumed.