29 - Sampling and Monte Carlo Methods: Estimating What You Cannot Compute

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 14

Goal. Turn simulation from a checking tool into a computing tool: estimate an expectation you cannot integrate, and know what the estimate is worth. By the end you can prove that the sample mean is unbiased and that its error falls as 1/√n, you can generate samples from any distribution by inverse transform or by rejection, you can move samples where they are needed with importance sampling, and you can explain why the reparametrisation trick makes differentiable an expectation that otherwise is not.

A myth to dispel. «Monte Carlo is an approximate method: if you want the right number, use a formula.» No: Monte Carlo is a method that is exact on average with a quantifiable error, and in high dimension it is the only one that works. The closed form, when it exists, always wins; the point is that it almost never exists.

Introduction

In the lesson on random variables simulation served as a check: numbers were generated to watch the law of large numbers and the central limit theorem «happen», that is, to convince ourselves of two theorems already proved. Here simulation changes role and becomes an estimator: the tool with which you compute an integral nobody can solve by hand.

This is not a side topic. Three later lessons use these tools as given: the variational autoencoder is trainable only thanks to reparametrisation, the instabilities of the ELBO are instabilities of a Monte Carlo estimator, and the whole reinforcement learning module, starting from the Bellman equation, estimates expected returns by simulation. Whoever arrives there without this lesson is copying formulas.

Final outcome of this lesson. Writing an estimator by simulation, declaring its standard error, choosing the number of samples from the error you want, and recognising the ways in which a correct estimator still produces unusable numbers.

Roadmap of the lesson.

1. The Monte Carlo estimator, the proof that it is unbiased, and its error σ/√n.
2. The comparison with grid quadrature, in dimension 1, 3 and 10.
3. How random numbers are generated, and why the seed is part of the experiment.
4. Inverse transform: derivation, continuous and discrete example, limits.
5. Rejection sampling: proof of correctness and collapse in high dimension.
6. Importance sampling, its variance, and the case where it blows up.
7. Antithetic variables and control variates, with the reduction measured.
8. The reparametrisation trick, compared numerically with REINFORCE.
9. A complete experiment on a five-state queue, with confidence intervals.

The outline goes from the method to its fragilities: first the estimator and its error, then the three ways of producing samples (inverse transform, rejection, importance), then how to reduce variance, finally reparametrisation and a complete experiment. Every technique has a section in which it breaks: read those before trusting an estimate.

The Monte Carlo idea, and the price of its precision

The analogy: counting the fish in the lake

To know how many fish there are in a lake you do not drain the lake. You catch a hundred fish, tag them and put them back; a few days later you fish again and look at what fraction of the catch carries a tag. If twelve out of a hundred recaught fish are tagged, the hundred tagged ones are about 12% of the population: eight hundred fish.

Two things about that estimate are worth noting, and they are exactly the two things of Monte Carlo. The method does not look at the whole lake, it looks at a sample, and it works precisely because the sample is taken at random. And the result is not a number, it is a number with an uncertainty: fishing again you would find eleven tagged, or fourteen. A biologist reporting «eight hundred fish» without saying «give or take how many» would not have made a measurement. In engineering terms: when a quantity is written as an expected value, you can estimate it by drawing samples and averaging them.

The estimator, and the proof that it is unbiased

Let X be a random variable with density p, and f a function whose mean value is of interest. The quantity you want is

I = E_p[f(X)] = ∫ f(x) p(x) dx.

The Monte Carlo estimator takes n independent samples x₁, …, xₙ drawn from p and averages them:

Î_n = (1/n) · Σᵢ f(xᵢ).

That this estimator is unbiased is proved in two lines, and the proof uses only the linearity of expectation: E[Î_n] = E[(1/n) Σᵢ f(xᵢ)] = (1/n) Σᵢ E[f(xᵢ)] = (1/n) · n · I = I. In words: the mean of the means is the mean. There is no approximation in this step, and no assumption on f beyond the existence of I. ∎

«Unbiased» means that on average the estimator gives the right answer, not that one run does. These are two different statements, and confusing them is mistake number one in this subject.

The error: σ/√n, and why one decimal costs a hundred times

How far is Î_n from I? The f(xᵢ) are independent and all have variance σ² = Var_p[f(X)]; the variance of a sum is the sum of the variances, and dividing by n squares the divisor:

Var[Î_n] = (1/n²) · Σᵢ Var[f(xᵢ)] = (1/n²) · n · σ² = σ²/n,

so the standard deviation of the estimator, which is called the standard error, is σ/√n. Read it as a sentence: the uncertainty falls as the square root of the number of samples. The consequence is brutal and must be said explicitly: to gain one decimal digit you need a hundred times the samples. From n to 4n the error halves, from n to 100n it is divided by ten. Whoever promises «just simulate a bit more» has not done this arithmetic.

There is, however, a property that makes Monte Carlo irreplaceable, and it lies in what σ/√n does not contain: the dimension of the space. Integrating in one dimension or in a thousand does not change the rate of convergence, it only changes σ. No deterministic method has this property, and that is why Bayesian statistics, computational physics and reinforcement learning all run on simulations.

Let us check both claims on a case with a known answer: X ∼ N(0, 1) and f(x) = x², so that I = 1 and σ² = Var[X²] = 2, which makes the theoretical standard error √2/√n.

import numpy as np

def stima_monte_carlo(n, rng):
    """Estimate E[f(X)] with f(x) = x² and X ~ N(0, 1): the exact value is 1."""
    x = rng.standard_normal(n)
    f = x ** 2
    return f.mean(), f.std(ddof=1) / np.sqrt(n)      # σ̂/√n, NON σ̂

rng = np.random.default_rng(42)
print("      n   estimate  estimated error    true error   σ/√n theory")
for n in (10, 100, 1000, 10000, 100000):
    stima, errore_stimato = stima_monte_carlo(n, rng)
    print(f"{n:>7} {stima:>10.5f} {errore_stimato:>16.5f}"
          f" {abs(stima - 1.0):>13.5f} {np.sqrt(2.0 / n):>13.5f}")

# the typical error is measured by repeating the experiment, not by running it once
rng = np.random.default_rng(7)
print("\n      n    root mean square error    σ/√n theory      ratio")
for n in (10, 100, 1000, 10000, 100000):
    errori = np.array([(rng.standard_normal(n) ** 2).mean() - 1.0 for _ in range(1000)])
    rmse = float(np.sqrt((errori ** 2).mean()))
    print(f"{n:>7} {rmse:>25.6f} {np.sqrt(2.0 / n):>14.6f} {rmse / np.sqrt(2.0 / n):>10.4f}")

Output obtained (Python 3.14.4, NumPy 2.3.5, seeds 42 and 7):

      n   estimate  estimated error    true error   σ/√n theory
     10    0.89688          0.36726       0.10312       0.44721
    100    0.59126          0.07694       0.40874       0.14142
   1000    1.00197          0.04563       0.00197       0.04472
  10000    1.01770          0.01452       0.01770       0.01414
 100000    1.00663          0.00452       0.00663       0.00447

      n    root mean square error    σ/√n theory      ratio
     10                  0.431696       0.447214     0.9653
    100                  0.140293       0.141421     0.9920
   1000                  0.044682       0.044721     0.9991
  10000                  0.014107       0.014142     0.9975
 100000                  0.004492       0.004472     1.0043

What the code does. The first table runs one simulation for each n: estimate, standard error from the data (σ̂/√n, where the /√n is the part almost everyone forgets), true error — which we know only because the answer is known — and theoretical error. The second repeats the experiment a thousand times for each n and reports the root mean squared error: this is how the typical error of an estimator is measured.

What it tells us. The second table confirms the law across four orders of magnitude: the ratio between measured error and theoretical σ/√n lies between 0.965 and 1.004. The first says something more interesting at row n = 100, where the true error (0.409) is more than five times the estimated one (0.077). It is not a bug: it is an unlucky sample in which the estimate is low (0.591) and the sample variance is underestimated along with it, because x² is skewed and whoever misses the right tail misses mean and spread together. A standard error is itself an estimate, on small n it is as noisy as the quantity it watches over, and it is moreover correlated with the error it is supposed to measure.

Monte Carlo versus a grid: when it really pays

The deterministic alternative to Monte Carlo is grid quadrature: you split the domain into cells, evaluate the function at the centre of each and sum. In one dimension the midpoint rule has error O(h²) with h = 1/m the step, that is O(m⁻²): doubling the nodes divides the error by four. Monte Carlo, by comparison, with its O(n⁻¹/²) cuts a very poor figure.

The trouble is that in d dimensions a grid with m nodes per axis needs m^d of them. For the same budget of evaluations N, the grid can afford only m ≈ N^{1/d} nodes per axis, so its error is O(N^{−2/d}). Compare it with O(N^{−1/2}): the grid wins while 2/d > 1/2, that is while d < 4, and loses afterwards. It is the curse of dimensionality seen in a one-line inequality.

Let us measure it on ∫ exp(−‖x‖²) dx over the cube [0,1]^d, which is separable and therefore has exact value ((√π/2)·erf(1))^d, with a fixed budget of 10 000 evaluations.

import math
import numpy as np

def integrale_esatto(d):
    """∫ exp(−‖x‖²) over [0,1]^d: the integrand is separable, so the result is a power."""
    return (math.sqrt(math.pi) / 2.0 * math.erf(1.0)) ** d

def quadratura_punto_medio(d, bilancio):
    """Midpoint rule: m nodes per axis, hence m^d evaluations."""
    m = int(bilancio ** (1.0 / d))
    nodi = (np.arange(m) + 0.5) / m
    una_dimensione = float(np.exp(-nodi ** 2).sum() / m)
    return una_dimensione ** d, m, m ** d

def quadratura_monte_carlo(d, n, rng):
    u = rng.random((n, d))
    f = np.exp(-(u ** 2).sum(axis=1))
    return f.mean(), f.std(ddof=1) / np.sqrt(n)

rng = np.random.default_rng(42)
bilancio = 10_000
print("  d | m per axis  evaluations      grid error |  MC error   MC standard error")
for d in (1, 3, 10):
    esatto = integrale_esatto(d)
    griglia, m, usate = quadratura_punto_medio(d, bilancio)
    mc, errore_standard = quadratura_monte_carlo(d, bilancio, rng)
    print(f" {d:>2} | {m:>10} {usate:>12} {abs(griglia - esatto):>15.3e} |"
          f" {abs(mc - esatto):>9.3e} {errore_standard:>19.3e}")

Output obtained (seed 42):

  d | m per axis  evaluations      grid error |  MC error   MC standard error
  1 |      10000        10000       3.066e-10 | 2.105e-03           2.005e-03
  3 |         21         9261       1.163e-04 | 3.780e-03           1.998e-03
 10 |          2         1024       5.889e-03 | 3.253e-05           5.457e-04

What the code does. quadratura_punto_medio exploits separability: the integral over d dimensions is the d-th power of the one-dimensional one. int(bilancio ** (1.0 / d)) is the number of nodes per axis the budget allows: 10 000 in one dimension, 21 in three (21³ = 9261), and only 2 in ten, because 3¹⁰ = 59 049 would already be six times over budget.

What it tells us, and it must be read honestly. In dimension 1 the grid wins outright: 3.07·10⁻¹⁰ against 2.11·10⁻³, seven orders of magnitude for the same evaluations. Anyone integrating a smooth function of one variable with Monte Carlo is wasting a factor of ten million. In dimension 3 the grid is still thirty times ahead. In dimension 10 the ratio flips: throttled to two nodes per axis, the grid is off by 5.9·10⁻³ against a Monte Carlo standard error of 5.5·10⁻⁴.

A warning about the d = 10 row. The Monte Carlo error reported there is 3.25·10⁻⁵, that is smaller than its own standard error 5.46·10⁻⁴: it is a lucky draw, the same phenomenon as the n = 100 row above with the sign reversed. The number to compare with the grid is 5.46·10⁻⁴; reporting the other would mean claiming a precision the method does not have. The rule to take away: Monte Carlo is not «better» than quadrature, it is indifferent to dimension. Below four dimensions, with smooth functions, use a grid; above, the grid does not exist as an option.

Generating random numbers: the seed is part of the experiment

Everything so far assumes you can draw samples, and at the bottom there is a single primitive: a generator of uniform numbers in [0, 1). All the rest — normals, exponentials, discrete distributions, permutations — is built from there, and that is what the next two sections do.

Those numbers are not random: they are produced by a deterministic function of the generator's internal state, whose initial value is the seed. They are called pseudorandom because they pass the statistical tests of randomness while being reproducible, and this is not a weakness to tolerate: it is the property that makes a Monte Carlo experiment an experiment rather than an anecdote. The seed is part of the measurement protocol and must be declared the way a library version is declared.

In NumPy you use np.random.default_rng(seme), which returns a generator with its own state, based on O'Neill's PCG64 algorithm. The old interface np.random.seed() plus np.random.rand() is declared frozen by the official documentation, which reserves it for cases where you need the same values as old versions, and it should be avoided in new code because it acts on a global state shared with every library in the process: it only takes one dependency drawing a number to shift your sequence, and the «reproducible» result stops being so with nothing to flag it.

import numpy as np

rng_a = np.random.default_rng(42)
rng_b = np.random.default_rng(42)
print("same seed, two generators:", rng_a.random(3), "|", rng_b.random(3))

rng = np.random.default_rng(42)
print("first block: ", rng.random(3))
print("second block:", rng.random(3))       # the state has advanced: different numbers
print("algorithm:", type(rng.bit_generator).__name__)

# two independent, reproducible streams, without touching any global state
radice = np.random.SeedSequence(42)
flussi = [np.random.default_rng(s) for s in radice.spawn(2)]
print("stream 0:", flussi[0].random(2))
print("stream 1:", flussi[1].random(2))

Output obtained (seed 42):

same seed, two generators: [0.77395605 0.43887844 0.85859792] | [0.77395605 0.43887844 0.85859792]
first block:  [0.77395605 0.43887844 0.85859792]
second block: [0.69736803 0.09417735 0.97562235]
algorithm: PCG64
stream 0: [0.91674416 0.91098667]
stream 1: [0.46749078 0.0464489 ]

What the code does. Two generators with the same seed produce the same sequence digit for digit, but the state advances: the second call to the same generator gives different numbers. The last block uses SeedSequence.spawn, the correct way to obtain independent streams from a single seed when simulating in parallel — putting default_rng(42) in every process would repeat the exact same work, with zero variance reduction and no error. The rule to take away: a generator is passed as an argument. A function that calls np.random.rand() inside itself is not testable, and two functions that do it are not composable.

Inverse transform sampling

The derivation

Given a random variable with cumulative distribution function (CDF) F(x) = P(X ≤ x) continuous and strictly increasing, you want to generate samples of X starting from U ∼ U(0,1). The claim is that applying F⁻¹ is enough:

if U ∼ U(0,1) then X = F⁻¹(U) has cumulative distribution function F.

The proof is three steps and uses only the definition. Compute the CDF of F⁻¹(U): P(F⁻¹(U) ≤ x). Since F is increasing, applying it to both sides does not change the direction of the inequality, so the event {F⁻¹(U) ≤ x} coincides with {U ≤ F(x)}. Finally, for a uniform on [0,1] we have P(U ≤ t) = t for every t ∈ [0,1], and F(x) lies right inside that range. Putting it together: P(F⁻¹(U) ≤ x) = P(U ≤ F(x)) = F(x), which is the claim. ∎

The geometric intuition: where the density is high the CDF rises steeply, so an interval of u of fixed width is mapped to a narrow interval of x, and many samples fall there. Inverse transform redistributes uniform points according to the slope of F.

Two examples worked through to the end

Exponential. For X ∼ Exp(λ) we have F(x) = 1 − e^{−λx} for x ≥ 0. Set u = 1 − e^{−λx} and solve for x: e^{−λx} = 1 − u, hence −λx = ln(1 − u), hence x = −ln(1 − u)/λ. One line of code, no rejection, no loop.

Discrete. With outcomes a₁, …, a_k and probabilities p₁, …, p_k the CDF is a staircase: compute the cumulatives c_j = p₁ + … + p_j and return the first a_j with u ≤ c_j. The generalised inverse is F⁻¹(u) = min{ a_j : F(a_j) ≥ u }, in NumPy a single call to searchsorted.

import numpy as np

def esponenziale_per_trasformata_inversa(n, lam, rng):
    """F⁻¹(u) = −ln(1 − u)/λ. `rng.random()` lies in [0, 1): 1 − u is never 0."""
    u = rng.random(n)
    return -np.log(1.0 - u) / lam

rng = np.random.default_rng(42)
lam = 0.4
campioni = esponenziale_per_trasformata_inversa(200_000, lam, rng)
print(f"mean     {campioni.mean():.5f}   expected 1/λ    = {1 / lam:.5f}")
print(f"median   {np.median(campioni):.5f}   expected ln2/λ  = {np.log(2) / lam:.5f}")
print(f"std.dev. {campioni.std(ddof=1):.5f}   expected 1/λ    = {1 / lam:.5f}")

# discrete case: the outcome of a service request, three possible values
esiti = ("served", "queued", "rejected")
probabilita = np.array([0.6, 0.3, 0.1])
cumulate = np.cumsum(probabilita)
indici = np.searchsorted(cumulate, rng.random(200_000))
for k, nome in enumerate(esiti):
    print(f"{nome:>10}: frequency {np.mean(indici == k):.5f}   expected {probabilita[k]:.5f}")

Output obtained (seed 42):

mean     2.50114   expected 1/λ    = 2.50000
median   1.73343   expected ln2/λ  = 1.73287
std.dev. 2.50291   expected 1/λ    = 2.50000
    served: frequency 0.60023   expected 0.60000
    queued: frequency 0.30019   expected 0.30000
  rejected: frequency 0.09959   expected 0.10000

What the code does. The function applies literally the formula derived above. The detail of 1.0 - u instead of u is not cosmetic: rng.random() gives values in [0, 1), so zero is possible and one is not, and 1 − u never vanishes; writing −ln(u)/λ would produce an inf with tiny but non-zero probability, that is a failure that shows up in production after weeks. In the discrete case searchsorted returns the index of the first cumulative that exceeds u: the generalised inverse in one line, in O(log k).

What it tells us. All three statistics of the exponential and the three discrete frequencies land on the third digit. That mean and standard deviation both equal 1/λ is a property of the exponential, not a typo.

The limit of the method, and it is a big one. You need F⁻¹ in closed form. For the exponential, the uniform, the Weibull, the Cauchy and finite discrete ones it exists; for the normal it does not (its CDF is a special function, invertible only numerically), and for almost every distribution of interest in Bayesian statistics it does not either.

Rejection sampling

The geometric idea

Sampling from a density p is the same as drawing a uniform point in the area under its graph and keeping the abscissa. Rejection sampling says: if I cannot draw points under p, I draw them under a more convenient curve that lies above it, and throw away those that land in the space between the two. You need a proposal q you can sample from and a constant M such that

p(x) ≤ M · q(x) for every x.

The algorithm is: draw x ∼ q; draw u ∼ U(0,1); accept x if u ≤ p(x) / (M·q(x)), otherwise start again.

The proof that accepted samples have distribution p

Let us compute the joint probability of drawing a value ≤ x and accepting it. Conditioning on the drawn value: P(X ≤ x, accetto) = ∫_{−∞}^{x} q(t) · [p(t)/(M q(t))] dt, because the probability of accepting given t is exactly that ratio. The two q(t) cancel and what remains is (1/M) ∫_{−∞}^{x} p(t) dt = F_p(x)/M.

Letting x → +∞ in the same formula gives the probability of accepting, whatever the value: P(accetto) = 1/M. Dividing the two gives the conditional distribution of the accepted samples alone: P(X ≤ x | accetto) = [F_p(x)/M] / (1/M) = F_p(x). ∎

Two operational consequences. The accepted samples are exact, whichever valid M is used; and the acceptance rate is 1/M, so on average M attempts are needed per sample. A needlessly large M introduces no bias, it only costs time — and that is a subtle nuisance, because a slow but correct method raises no exceptions.

Let us check it on a Beta(2, 5), whose density 30·x·(1−x)⁴ has no closed-form inverse CDF, with a uniform proposal. The maximum is at x = 1/5 and equals M = 2.4576, so the theoretical rate is 1/M = 0.406901.

import numpy as np
from scipy.stats import beta

def densita_bersaglio(x):
    """Beta(2, 5) density written by hand: 30·x·(1 − x)⁴ on [0, 1]."""
    return 30.0 * x * (1.0 - x) ** 4

M = densita_bersaglio(0.2)        # maximum of the density, with uniform proposal q(x) = 1

def campiona_per_rifiuto(tentativi, rng):
    x = rng.random(tentativi)                     # proposal q = U(0, 1)
    altezza = rng.random(tentativi) * M           # uniform point under M·q(x)
    accettati = altezza <= densita_bersaglio(x)
    return x[accettati], accettati.mean()

rng = np.random.default_rng(42)
campioni, tasso = campiona_per_rifiuto(500_000, rng)
print(f"M = {M:.6f}   theoretical acceptance rate 1/M = {1 / M:.6f}")
print(f"measured rate {tasso:.6f}   accepted samples {campioni.size}")
print(f"mean     {campioni.mean():.6f}   expected 2/7 = {2 / 7:.6f}")
print(f"variance {campioni.var(ddof=1):.6f}   expected     = {beta.var(2, 5):.6f}")
for q in (0.1, 0.5, 0.9):
    print(f"quantile {q:.1f}: measured {np.quantile(campioni, q):.6f}"
          f"   exact {beta.ppf(q, 2, 5):.6f}")

Output obtained (seed 42, SciPy 1.16.3):

M = 2.457600   theoretical acceptance rate 1/M = 0.406901
measured rate 0.407112   accepted samples 203556
mean     0.285713   expected 2/7 = 0.285714
variance 0.025444   expected     = 0.025510
quantile 0.1: measured 0.093102   exact 0.092595
quantile 0.5: measured 0.264270   exact 0.264450
quantile 0.9: measured 0.509763   exact 0.510316

What the code does. Instead of drawing u ∼ U(0,1) and comparing it with p(x)/(M·q(x)), it draws a uniform height in [0, M] and compares it with p(x): it is the same inequality multiplied by M, but it reads for what it is — a uniform point in the rectangle, whose abscissa is kept if it falls under the curve. The vectorised version generates all attempts at once and masks, rather than using a while: you pay in memory what you gain in speed, and the number of samples obtained is not known in advance.

What it tells us. Measured rate 0.407112 against a theoretical 0.406901, five digits of agreement, and the 203 556 samples accepted out of 500 000 are exactly what the rate predicts. Mean, variance and three quantiles match SciPy's exact values to the third or fourth digit.

Why rejection stops working in high dimension

The acceptance rate 1/M is a ratio of volumes, and ratios of volumes in high dimension behave counterintuitively. The sharpest experiment is the simplest: draw uniform points in the cube [−1, 1]^d and accept those inside the unit ball. The rate is the ratio between the volume of the ball, π^{d/2}/Γ(d/2 + 1), and that of the cube, 2^d.

import numpy as np
from scipy.special import gammaln

def tasso_teorico(d):
    """Volume of the unit ball divided by the volume of the cube [−1, 1]^d."""
    log_palla = (d / 2) * np.log(np.pi) - gammaln(d / 2 + 1)
    return float(np.exp(log_palla - d * np.log(2.0)))

rng = np.random.default_rng(42)
n = 200_000
print("  d |  measured rate     theory rate       tries per 1 sample")
for d in (1, 2, 5, 10, 20):
    punti = rng.uniform(-1.0, 1.0, size=(n, d))
    dentro = (punti ** 2).sum(axis=1) <= 1.0
    print(f" {d:>2} | {dentro.mean():>14.6f} {tasso_teorico(d):>15.8f}"
          f" {1 / tasso_teorico(d):>24.1f}")

Output obtained (seed 42):

  d |  measured rate     theory rate       tries per 1 sample
  1 |       1.000000      1.00000000                      1.0
  2 |       0.786720      0.78539816                      1.3
  5 |       0.165875      0.16449341                      6.1
 10 |       0.002700      0.00249039                    401.5
 20 |       0.000000      0.00000002               40631627.6

What the code does. gammaln computes the logarithm of the Gamma function: you work on a logarithmic scale because π^{d/2} and Γ(d/2+1) both become enormous while their ratio becomes tiny, and in double precision the numerator overflows long before the result loses meaning. The comparison is on the squared norm, which avoids a useless square root.

What it tells us. You accept 79% of the attempts in two dimensions, 16% in five, 0.25% in ten — already 400 attempts per sample. In twenty the theoretical rate is 2.5·10⁻⁸ and the simulation, over 200 000 attempts, accepts zero: on average 40 million would be needed for a single sample. No optimisation fixes this: it is the volume of the cube concentrating in the corners. The rule to take away: rejection is exact, simple and excellent in one or two dimensions; beyond five it is dead, and that is why Markov chain Monte Carlo exists, which does not draw independent samples but builds a path that visits the high-density regions.

Importance sampling

The rewriting, and the weight

Back in the lake of the opening analogy: if the fish you care about are all in one inlet, fishing at random over the whole lake is a waste, and it is better to fish where they are and then correct the count by how much you oversampled them. That is importance sampling, and the rewriting that justifies it fits on one line: you multiply and divide by a proposal density q of your choosing, provided it is positive wherever f·p is not zero.

E_p[f(X)] = ∫ f(x) p(x) dx = ∫ f(x) · [p(x)/q(x)] · q(x) dx = E_q[ f(X) · w(X) ], with w(x) = p(x)/q(x).

In words: you can estimate a mean with respect to p by drawing from any other distribution q, provided every contribution is corrected with the importance weight w = p/q, large where q undersamples and small where it oversamples. The positivity condition is not a quibble: if q is zero where f·p is not, that region is never visited and the estimator becomes biased with no symptom at all.

The variance is Var_q[f·w] = ∫ f(x)² p(x)²/q(x) dx − I², and that integral is the recipe for choosing q: put mass where |f|·p is large. The optimal proposal q*(x) ∝ |f(x)|·p(x) zeroes the variance for f ≥ 0, and is unusable because its normalising constant is the very integral you are after; it still serves as a compass.

The gain is spectacular on rare events. Let us estimate P(X > 4) for X ∼ N(0,1), which equals 3.167·10⁻⁵: with naive Monte Carlo, out of a hundred thousand samples, three are expected above the threshold. With proposal q = N(4, 1) half the samples exceed it, and the weight is p(x)/q(x) = exp(8 − 4x).

import numpy as np
from scipy.stats import norm

soglia = 4.0
esatto = float(norm.sf(soglia))          # P(X > 4) with X ~ N(0, 1)

def monte_carlo_ingenuo(n, rng):
    x = rng.standard_normal(n)
    indicatore = (x > soglia).astype(float)
    return indicatore.mean(), indicatore.std(ddof=1) / np.sqrt(n)

def monte_carlo_per_importanza(n, rng):
    """Proposal q = N(4, 1): the weight p/q is exp(8 − 4x)."""
    x = rng.standard_normal(n) + soglia
    contributo = np.exp(8.0 - 4.0 * x) * (x > soglia)
    return contributo.mean(), contributo.std(ddof=1) / np.sqrt(n)

rng = np.random.default_rng(42)
n = 100_000
print(f"exact value Φ(−4) = {esatto:.6e}")
for nome, stimatore in (("naive", monte_carlo_ingenuo),
                        ("importance", monte_carlo_per_importanza)):
    stima, errore_standard = stimatore(n, rng)
    relativo = errore_standard / stima if stima > 0 else float("inf")
    print(f"{nome:>15}: estimate {stima:.6e}   standard error {errore_standard:.3e}"
          f"   relative error {relativo:>8.2%}")

Output obtained (seed 42):

exact value Φ(−4) = 3.167124e-05
          naive: estimate 9.000000e-05   standard error 3.000e-05   relative error   33.33%
     importance: estimate 3.173136e-05   standard error 2.127e-07   relative error    0.67%

What the code does. The naive estimator counts how many samples exceed the threshold: it is the frequency of a rare event, that is a binomial with very few successes. The importance one shifts the proposal onto the threshold (rng.standard_normal(n) + soglia is a sample from N(4,1)) and weights every contribution by exp(8 − 4x), obtained by simplifying the ratio of the two Gaussian densities.

What it tells us. A relative error of 0.67% against 33%, that is a variance fifty times smaller for the same number of samples. And there is more: the naive estimate is 9.0·10⁻⁵, almost three times the true value, because it happened to catch nine samples beyond the threshold instead of the three expected; its 95% interval [3.0·10⁻⁵, 1.5·10⁻⁴] contains the true value by a hair. It is a correct estimator that, on a rare event, produces a number nobody should decide on.

When the variance blows up: tails that are too light

The way importance sampling fails is particularly unpleasant: the estimator stays unbiased but acquires infinite variance, and an infinite variance is not visible by looking at the estimate. The condition lies in ∫ f²p²/q: if q has lighter tails than f·p, the ratio p/q diverges in the tails faster than q damps them, and the integral does not converge.

In numbers: target p = N(0,1), f(x) = x², exact value 1, proposal q = N(0, σ_q²) with weight σ_q · exp(x²(1/(2σ_q²) − 1/2)). For σ_q = 1.5 the exponent is negative and all is well; for σ_q = 0.5 it is +1.5x² and the second moment contains ∫ x⁴ e^{x²} dx, which diverges.

import numpy as np

def contributi_per_importanza(n, sigma_proposta, rng):
    """Target p = N(0, 1) and f(x) = x², so E_p[f] = 1 exactly."""
    x = rng.standard_normal(n) * sigma_proposta
    log_peso = -x ** 2 / 2.0 + x ** 2 / (2.0 * sigma_proposta ** 2) + np.log(sigma_proposta)
    return x ** 2 * np.exp(log_peso)

rng = np.random.default_rng(42)
for sigma in (1.5, 0.5):
    contributi = contributi_per_importanza(200_000, sigma, rng)
    corse = np.cumsum(contributi) / np.arange(1, contributi.size + 1)
    print(f"proposal N(0, {sigma}²)   tails {'heavier' if sigma > 1 else 'lighter'} than the target")
    for n in (1_000, 10_000, 50_000, 100_000, 200_000):
        print(f"   n = {n:>7}: estimate {corse[n - 1]:.6f}")
    print(f"   largest contribution observed {contributi.max():.3e}"
          f"   sample variance {contributi.var(ddof=1):.3e}")

Output obtained (seed 42):

proposal N(0, 1.5²)   tails heavier than the target
   n =    1000: estimate 0.967729
   n =   10000: estimate 0.990041
   n =   50000: estimate 0.996901
   n =  100000: estimate 0.997347
   n =  200000: estimate 0.998902
   largest contribution observed 1.987e+00   sample variance 4.913e-01
proposal N(0, 0.5²)   tails lighter than the target
   n =    1000: estimate 0.904171
   n =   10000: estimate 0.921567
   n =   50000: estimate 0.812826
   n =  100000: estimate 0.903204
   n =  200000: estimate 0.847854
   largest contribution observed 9.829e+03   sample variance 6.363e+02

What the code does. The weight is computed on a logarithmic scale and exponentiated once only: with σ_q = 0.5 the exponent reaches 1.5·x², which for x ≈ 2.3 is already e⁸. np.cumsum divided by the indices produces the running mean: the convergence plot, read as a table.

What it tells us, and this is the point of the section. With the wide proposal the running mean climbs tidily towards 1 (0.9677 → 0.9900 → 0.9969 → 0.9973 → 0.9989), maximum contribution 1.99 and variance 0.49. With the narrow proposal it does not converge: 0.904 → 0.922 → 0.813 → 0.903 → 0.848. Look at the step from 10 000 to 50 000: the estimate gets worse on five times the samples. The culprit is the maximum contribution, 9.8·10³, five thousand times that of the good proposal: a single sample is worth more than all the others together, and the mean jumps every time one arrives.

The symptom and the remedy. An estimate that seems to converge and then jumps, or a sample variance that grows with n, almost always means proposal tails lighter than the target. Always print the maximum weight and the effective sample size ESS = (Σwᵢ)²/Σwᵢ²: if over two hundred thousand samples the ESS is in the tens, you do not have two hundred thousand samples, you have a few dozen. In case of doubt choose a proposal with heavier tails than the target: a weight tending to zero in the tails costs efficiency, a weight that blows up costs the validity of the estimate.

Variance reduction: antithetic variables and control variates

Since the error falls as 1/√n, increasing n is the most expensive way to become precise. The alternative is to lower σ: build an estimator of the same quantity with smaller variance. Two techniques cost almost nothing.

Antithetic variables. Instead of n independent samples you use n/2 pairs (U, 1−U). The variance of a pair mean is Var[(f(U) + f(1−U))/2] = [Var f + Cov]/2, because the two terms have the same distribution. If f is monotone, f(U) and f(1−U) are negatively correlated, Cov < 0, and the variance drops — for free, because the pair costs one random number less, not more.

Control variates. You exploit a g whose mean is known exactly, estimating f(x) − c·(g(x) − E[g]). The added term has zero mean, so the estimator stays unbiased for any c, which acts only on the variance Var f − 2c·Cov(f,g) + c²·Var g. That is a parabola in c, minimal at c* = Cov(f,g)/Var(g); substituting gives Var f · (1 − ρ²). The more g resembles f, the more you gain.

Let us measure them on f(u) = e^u with u uniform in [0,1] — the cost of an operation growing exponentially with load, synthetic data invented by us — whose exact value is e − 1. As control variate, g(u) = u with E[g] = 1/2. The comparison is at equal evaluations of f: otherwise the antithetic one would cheat.

import math
import numpy as np

def costo(u):
    """Cost of an operation that grows exponentially with the load u ∈ [0, 1]."""
    return np.exp(u)

esatto = math.e - 1.0
rng = np.random.default_rng(42)
n = 200_000                           # evaluations of f: the same for the three estimators

semplice = costo(rng.random(n))

v = rng.random(n // 2)
antitetico = (costo(v) + costo(1.0 - v)) / 2.0

u = rng.random(n)
c = float(np.cov(costo(u), u, ddof=1)[0, 1] / u.var(ddof=1))
controllo = costo(u) - c * (u - 0.5)

print(f"exact value e − 1 = {esatto:.6f}")
print("           estimator    estimate    variance per evaluation    standard error")
teoriche = {"plain": 0.2420356, "antithetic": 0.0078250, "control variate": 0.0039401}
for nome, campione, valutazioni in (("plain", semplice, 1),
                                    ("antithetic", antitetico, 2),
                                    ("control variate", controllo, 1)):
    varianza = campione.var(ddof=1) * valutazioni
    print(f"{nome:>20} {campione.mean():>11.6f} {varianza:>26.6f}"
          f" {campione.std(ddof=1) / np.sqrt(campione.size):>17.2e}"
          f"   (theoretical {teoriche[nome]:.6f})")
print(f"estimated coefficient c = {c:.6f}   theoretical optimum 12·Cov(e^U, U) = {12 * (1 - 0.5 * esatto):.6f}")

Output obtained (seed 42):

exact value e − 1 = 1.718282
           estimator    estimate    variance per evaluation    standard error
               plain    1.718252                   0.242112          1.10e-03   (theoretical 0.242036)
          antithetic    1.718174                   0.007763          1.97e-04   (theoretical 0.007825)
     control variate    1.718315                   0.003945          1.40e-04   (theoretical 0.003940)
estimated coefficient c = 1.689440   theoretical optimum 12·Cov(e^U, U) = 1.690309

What the code does. The three variances are reported per evaluation of f: the antithetic estimator uses two of them for every pair mean, so its variance must be multiplied by two before the comparison, otherwise the comparison is rigged in its favour by a factor of exactly 2. The coefficient c is estimated from the data with np.cov, as is done when Cov(f,g) cannot be computed by hand.

What it tells us. At equal 200 000 evaluations: 0.2421 for the plain estimator, 0.00776 for the antithetic one (31 times less), 0.00395 for the control variate (61 times less). Translated: to match by brute force the precision the control variate reaches with 200 000 evaluations you would need twelve million. The three measured variances hit the theoretical ones on the third digit, and the estimated coefficient (1.68944) hits the optimum (1.690309).

Trade-off. Antithetic variables require f to be monotone, or at least mostly so: on an f symmetric about 1/2 the correlation becomes positive and the variance increases. Control variates require a correlated g of known mean, which does not always exist; and estimating c on the same samples used for the mean introduces a bias that is negligible on large n but not on small n, where a separate pilot sample is preferable.

The reparametrisation trick

The problem: differentiating with respect to whoever generates the samples

So far the distribution being sampled from was given. In almost all modern machine learning it is not: you want to optimise it. That is, you need the gradient of

L(θ) = E_{x ∼ p_θ}[ f(x) ] with respect to θ.

And here the swap everybody makes — derivative inside the expectation — is not legitimate. ∇_θ E_{x∼p_θ}[f(x)] is not E_{x∼p_θ}[∇_θ f(x)], and the reason shows up when you write ∇_θ ∫ f(x) p_θ(x) dx: the dependence on θ is not in f, which often does not contain θ at all, it is in the measure. Differentiating f differentiates the wrong thing, and gives zero.

First solution: the score function estimator (REINFORCE)

You differentiate under the integral sign with respect to the density: ∇_θ ∫ f p_θ dx = ∫ f ∇_θ p_θ dx. But ∇_θ p_θ is not an expectation, because it is not something · p_θ. The fix is the identity ∇_θ p_θ = p_θ · ∇_θ log p_θ, which is the derivative of the logarithm read backwards:

∇_θ L(θ) = E_{x ∼ p_θ}[ f(x) · ∇_θ log p_θ(x) ].

This is the score function estimator, known in reinforcement learning as REINFORCE (Williams, 1992). It always works, even with f non-differentiable or discrete, because f is never differentiated. The price is heavy: f(x) multiplies a zero-mean noisy term, so the variance is high and grows with the absolute value of f. It is mitigated by subtracting a constant b from f (the baseline) — which is exactly a control variate — but it stays high.

Second solution: reparametrisation

The idea is to move the randomness out of the gradient path. If x ∼ N(μ, σ²) you can write x = μ + σ·ε with ε ∼ N(0,1): the distribution of ε does not depend on θ = (μ, σ), and all the dependence on the parameter sits in a deterministic function. Then

∇_θ E_{x∼p_θ}[f(x)] = ∇_θ E_{ε∼N(0,1)}[ f(μ + σε) ] = E_{ε}[ ∇_θ f(μ + σε) ],

and now the swap is legitimate, because the measure being integrated over no longer depends on θ: for μ you get E_ε[f′(μ + σε)], for σ you get E_ε[ε·f′(μ + σε)]. This is the reparametrisation trick (Kingma and Welling, 2014), and it is what makes a variational autoencoder trainable: without it, sampling in the latent space would break backpropagation. The price is that you need f differentiable and an explicit reparametrisation — conditions the Gaussian satisfies and a discrete distribution does not.

Let us compare the variances of the two estimators on ∇_μ E_{x∼N(μ,1)}[x²] with μ = 1.5, whose true value is 2μ = 3. The theoretical variances are computed by hand: μ⁴ + 14μ² + 15 = 51.5625 for the score function, 28 for the score with baseline, 4σ² = 4 for reparametrisation.

import numpy as np

mu, sigma = 1.5, 1.0
gradiente_vero = 2.0 * mu                # d/dμ E[x²] = d/dμ (μ² + σ²) = 2μ

def stimatore_score_function(n, rng):
    """REINFORCE: f(x) · ∂log p(x)/∂μ, with ∂log p/∂μ = (x − μ)/σ²."""
    x = mu + sigma * rng.standard_normal(n)
    return x ** 2 * (x - mu) / sigma ** 2

def stimatore_score_con_baseline(n, rng):
    """The same, with E[f] subtracted from f: it is a control variate."""
    x = mu + sigma * rng.standard_normal(n)
    return (x ** 2 - (mu ** 2 + sigma ** 2)) * (x - mu) / sigma ** 2

def stimatore_riparametrizzato(n, rng):
    """x = μ + σε with ε ∼ N(0,1): ∂f/∂μ = f′(x) · ∂x/∂μ = 2(μ + σε)."""
    eps = rng.standard_normal(n)
    return 2.0 * (mu + sigma * eps)

rng = np.random.default_rng(42)
teoriche = {"score function": mu ** 4 + 14 * mu ** 2 + 15,
            "score + baseline": 28.0,
            "reparametrisation": 4.0 * sigma ** 2}
print(f"true gradient = {gradiente_vero:.6f}")
print("           estimator   mean over 2·10⁶   measured variance    theory variance")
for nome, f in (("score function", stimatore_score_function),
                ("score + baseline", stimatore_score_con_baseline),
                ("reparametrisation", stimatore_riparametrizzato)):
    campione = f(2_000_000, rng)
    print(f"{nome:>20} {campione.mean():>17.5f} {campione.var(ddof=1):>19.4f}"
          f" {teoriche[nome]:>18.4f}")

print("\nstandard error of the gradient, at equal n")
print("      n   score function     reparametrisation")
for n in (10, 100, 1000, 10000):
    a = stimatore_score_function(n, rng)
    b = stimatore_riparametrizzato(n, rng)
    print(f"{n:>7} {a.std(ddof=1) / np.sqrt(n):>16.5f}"
          f" {b.std(ddof=1) / np.sqrt(n):>21.5f}")

Output obtained (seed 42):

true gradient = 3.000000
           estimator   mean over 2·10⁶   measured variance    theory variance
      score function           3.00276             51.7719            51.5625
    score + baseline           2.99788             27.9164            28.0000
   reparametrisation           3.00128              3.9978             4.0000

standard error of the gradient, at equal n
      n   score function     reparametrisation
     10          1.57648               0.51661
    100          1.02247               0.20037
   1000          0.25905               0.06371
  10000          0.07301               0.01981

What the code does. No automatic differentiation library: the two gradients are written by hand, and that is the best way to see that they are two different formulas for the same number. The score function multiplies f(x) = x² by the score (x − μ)/σ²; reparametrisation differentiates f and evaluates it at μ + σε. Both estimate 3, and they do: 3.00276 and 3.00128.

What it tells us. The measured variances hit the theoretical ones on the third digit (51.77 against 51.5625; 3.998 against 4.000), which confirms derivation and code together. The ratio is 12.9 in favour of reparametrisation, and the baseline recovers less than half the gap (51.6 → 28.0): useful, not a substitute. The second table translates variance into what you feel while training: with ten samples the score function has a standard error of 1.58 on a true value of 3, that is 53% — the sign is right, the magnitude is an opinion; reparametrisation sits at 0.52, that is 17%.

Trade-off, and it is sharp. Reparametrisation has far smaller variance but requires f differentiable and a reparametrisable distribution. The score function requires none of that, which is why it is the only usable one on discrete actions. It is not that one is better: they cover different cases, and the second is paid for with many more samples.

The experiment: a five-state queue

Let us put everything together on a small but non-trivial system, entirely synthetic: a service desk with a four-seat waiting room, hence five possible states (0, 1, 2, 3, 4 users in the system). Time is discrete; within each step, first the service and then the arrival: if there is at least one user, with probability s = 0.5 they are served and leave; after that, with probability a = 0.1 a new user arrives, who enters if there is room and is turned away if the system is full. The two probabilities are invented by us and declared as such: they exist to make the phenomena appear, not to describe a real desk.

The advantage of such a small system is that it has an exact answer to compare against, and a real system does not. The chain is birth-and-death: from a state n ≥ 1 you go up with probability (1−s)·a = 0.05 and down with s·(1−a) = 0.45, while from 0 you go up with a = 0.1. Detailed balance π_n · (up) = π_{n+1} · (down) gives the ratios π₁/π₀ = 0.1/0.45 and π_{n+1}/π_n = 0.05/0.45 = 1/9 for n ≥ 1; normalising yields the stationary distribution, and from it the mean number in the system L = Σ n·π_n. Little's law L = λ·W finally returns the mean wait W, where λ is the rate of accepted arrivals, that is a minus those turned away.

import numpy as np
from collections import deque

ARRIVO, SERVIZIO, CAPIENZA = 0.1, 0.5, 4

def riferimento_esatto():
    """Stationary distribution of the 5-state chain, and Little's law."""
    su = (1 - SERVIZIO) * ARRIVO                   # from n to n+1, for n ≥ 1
    giu = SERVIZIO * (1 - ARRIVO)                  # from n to n−1, for n ≥ 1
    rapporti = [1.0, ARRIVO / giu]                 # the first jump starts from 0
    for _ in range(3):
        rapporti.append(rapporti[-1] * su / giu)
    pi = np.array(rapporti) / sum(rapporti)
    numero_medio = float(np.arange(5) @ pi)
    arrivi_accettati = ARRIVO - pi[4] * su         # rejected arrivals do not enter
    return pi, numero_medio, numero_medio / arrivi_accettati

pi, numero_medio, attesa_esatta = riferimento_esatto()
print("stationary distribution:", np.round(pi, 6))
print(f"mean number in the system L = {numero_medio:.6f}"
      f"   exact mean wait W = L/λ = {attesa_esatta:.6f}")

Output obtained:

stationary distribution: [8.00024e-01 1.77783e-01 1.97540e-02 2.19500e-03 2.44000e-04]
mean number in the system L = 0.224851   exact mean wait W = L/λ = 2.248780

What the code does. It builds the ratios between consecutive stationary probabilities instead of solving the linear system: for a birth-and-death chain detailed balance gives them one at a time, and the final normalisation is a division by the sum. np.arange(5) @ pi is the dot product between state values and their probabilities, that is the expected value.

What it tells us. The desk is empty 80% of the time and full 0.024% of it. The exact mean wait is 2.2488 steps, to be compared with the 2 steps that service alone would cost if there were never a queue (the service duration is geometric with mean 1/s = 2): the queue adds barely a quarter of a step. This is the number the simulation must recover, and failing to recover it means the simulation is wrong.

Now the estimate by simulation, with the 95% confidence interval of the estimator — that is estimate ± 1.96 · standard error, where the standard error is computed on the spread of the per-repetition means, not on that of the individual waits. Two configurations: short runs of 400 steps, and long runs of 5000 steps of which the first 200 are discarded as warm-up.

def simula_coda(passi, rng, scarto=0):
    """Queue with 5 states (0..4): within each step first the service, then the arrival.
    A customer's wait is the number of step ends at which they are in the system.
    Only customers who arrived after `scarto` and left before the end are counted:
    those still inside at the last step are truncated, and must be discarded."""
    in_sistema, attese = deque(), []
    for t in range(passi):
        if in_sistema and rng.random() < SERVIZIO:
            arrivato = in_sistema.popleft()
            if arrivato >= scarto:
                attese.append(t - arrivato)
        if rng.random() < ARRIVO and len(in_sistema) < CAPIENZA:
            in_sistema.append(t)
    return np.array(attese, dtype=float)

rng = np.random.default_rng(42)
for etichetta, passi, scarto, lista in (
        ("short runs, no warm-up (400 steps)", 400, 0, (25, 100, 400, 1600)),
        ("long runs, 200 steps discarded (5000 steps)", 5000, 200, (25, 50, 100, 400))):
    print(f"\n{etichetta}")
    print("repetitions      mean wait    standard error   95% CI              covers W?")
    for ripetizioni in lista:
        medie = np.array([simula_coda(passi, rng, scarto).mean() for _ in range(ripetizioni)])
        stima = medie.mean()
        errore = medie.std(ddof=1) / np.sqrt(ripetizioni)
        basso, alto = stima - 1.96 * errore, stima + 1.96 * errore
        print(f"{ripetizioni:>11} {stima:>14.5f} {errore:>17.5f}"
              f"   [{basso:.5f}, {alto:.5f}]  {'yes' if basso <= attesa_esatta <= alto else 'NO'}")

Output obtained (seed 42):

short runs, no warm-up (400 steps)
repetitions      mean wait    standard error   95% CI              covers W?
         25        2.24170           0.06869   [2.10706, 2.37634]  yes
        100        2.25600           0.04361   [2.17052, 2.34147]  yes
        400        2.22397           0.01757   [2.18952, 2.25841]  yes
       1600        2.22711           0.00893   [2.20961, 2.24461]  NO

long runs, 200 steps discarded (5000 steps)
repetitions      mean wait    standard error   95% CI              covers W?
         25        2.24896           0.01883   [2.21206, 2.28586]  yes
         50        2.24007           0.01480   [2.21106, 2.26908]  yes
        100        2.23429           0.01117   [2.21240, 2.25617]  yes
        400        2.24940           0.00542   [2.23878, 2.26001]  yes

What the code does. The deque holds the arrival instants of the customers present, in order: service is FIFO, so the head of the queue always leaves first via popleft. The recorded wait is t − arrivato, that is the number of end-of-step counts in which that customer is in the system — exactly the quantity that appears in Little's law. The scarto parameter excludes from the count the customers who arrived during the warm-up. Customers still queueing at the last step are never recorded: it is an unavoidable truncation, and its consequence shows up in the first table.

What it tells us, and this is the central moment of the lesson. Look at the «copre W?» column of the first table. At 25, 100 and 400 repetitions the confidence interval contains the true value 2.24878. At 1600 repetitions it stops containing it: the interval is [2.20961, 2.24461] and the true value falls just outside. It is not bad luck and it is not an arithmetic slip: the 400-step runs produce an estimate biased downwards by about 0.02, for two reasons acting in the same direction — you start from an empty system and the first customers do not wait (transient), and the customers still inside at the end are not counted, although they are precisely the ones with the longest waits (truncation).

Read it carefully, because it is counterintuitive. While the repetitions were few, the interval was wide and the bias sat comfortably inside it: the protocol looked correct. As repetitions grew the interval tightened around the wrong number, and the bias, which had always been there, became visible only when precision grew enough to expose it. More data does not correct a bias: it only makes it more evident. The second table shows the remedy — runs twelve times longer with the first 200 steps thrown away — and there the interval covers the true value at all four sizes, including the tightest one (±0.0054).

Saturation: a rare event, and two ways of measuring it

A second question about the same system, this time about a rare event: starting from an empty system, what is the probability that the waiting room fills up (state 4) within eight steps? That needs at least four arrivals in eight steps, each with probability 0.1, and almost no service in between. Here too the exact answer exists, by dynamic programming on the chain with state 4 made absorbing.

import math
import numpy as np

ORIZZONTE = 8

def nucleo(peso_arrivo, peso_non_arrivo):
    """4×4 matrix on states 0..3 and the weight of the escape towards state 4.
    With (ARRIVO, 1−ARRIVO) it gives the true chain; with (a²/a′, (1−a)²/(1−a′))
    it gives the sum of p²/q, which is what the second moment needs."""
    Q = np.zeros((4, 4))
    Q[0, 0], Q[0, 1] = peso_non_arrivo, peso_arrivo
    for n in (1, 2, 3):
        Q[n, n - 1] = SERVIZIO * peso_non_arrivo
        Q[n, n] = SERVIZIO * peso_arrivo + (1 - SERVIZIO) * peso_non_arrivo
        if n + 1 <= 3:
            Q[n, n + 1] = (1 - SERVIZIO) * peso_arrivo
    return Q, (1 - SERVIZIO) * peso_arrivo

def somma_sui_percorsi(peso_arrivo, peso_non_arrivo):
    """Sum over the paths that touch state 4 within ORIZZONTE steps.
    `coda` is the factor for the steps that remain AFTER absorption: the weight
    keeps accumulating even when the event has already been decided."""
    Q, fuga = nucleo(peso_arrivo, peso_non_arrivo)
    coda = peso_arrivo + peso_non_arrivo
    u, totale = np.array([1.0, 0.0, 0.0, 0.0]), 0.0
    for k in range(ORIZZONTE):
        totale += fuga * u[3] * coda ** (ORIZZONTE - k - 1)
        u = u @ Q
    return totale

def varianza_teorica(arrivo_proposto):
    """Var_q[w·1_A] = E_q[w²·1_A] − P(A)², with E_q[w²·1_A] = Σ p²/q over the paths."""
    secondo = somma_sui_percorsi(ARRIVO ** 2 / arrivo_proposto,
                                 (1 - ARRIVO) ** 2 / (1 - arrivo_proposto))
    return secondo - esatto ** 2

def episodi(n, arrivo_proposto, rng):
    """n episodes of ORIZZONTE steps from an empty system, in parallel.
    Returns the contributions 1{touched state 4} · weight of the path."""
    stato = np.zeros(n, dtype=np.int64)
    saturato, log_peso = np.zeros(n, dtype=bool), np.zeros(n)
    sopra = math.log(ARRIVO / arrivo_proposto)
    sotto = math.log((1 - ARRIVO) / (1 - arrivo_proposto))
    for _ in range(ORIZZONTE):
        stato -= (stato > 0) & (rng.random(n) < SERVIZIO)
        arrivo = rng.random(n) < arrivo_proposto
        log_peso += np.where(arrivo, sopra, sotto)
        entra = arrivo & (stato < CAPIENZA)
        stato += entra
        saturato |= entra & (stato == CAPIENZA)
    return saturato * np.exp(log_peso)

esatto = somma_sui_percorsi(ARRIVO, 1 - ARRIVO)      # true weights: it is the probability
rng = np.random.default_rng(42)
print(f"exact probability of saturation in {ORIZZONTE} steps = {esatto:.6e}")
print("      estimator           n     estimate    standard error"
      "    relative error useful episodes    theory variance")
misurate = {}
for nome, n, arrivo_proposto in (("naive", 20_000, ARRIVO),
                                 ("naive", 2_000_000, ARRIVO),
                                 ("importance", 20_000, 0.5)):
    contributi = episodi(n, arrivo_proposto, rng)
    stima, varianza = contributi.mean(), contributi.var(ddof=1)
    misurate[nome] = varianza
    print(f"{nome:>15} {n:>11} {stima:>12.6e} {math.sqrt(varianza / n):>17.3e}"
          f" {math.sqrt(varianza / n) / stima if stima > 0 else float('nan'):>17.1%}"
          f" {int((contributi > 0).sum()):>15} {varianza_teorica(arrivo_proposto):>18.6e}")
print(f"variance ratio, measured    = "
      f"{misurate['naive'] / misurate['importance']:.1f}")
print(f"variance ratio, theoretical = "
      f"{varianza_teorica(ARRIVO) / varianza_teorica(0.5):.1f}")

Output obtained (seed 42):

exact probability of saturation in 8 steps = 2.217811e-04
      estimator           n     estimate    standard error    relative error useful episodes    theory variance
          naive       20000 2.000000e-04         9.999e-05             50.0%               4       2.217319e-04
          naive     2000000 2.185000e-04         1.045e-05              4.8%             437       2.217319e-04
     importance       20000 2.125971e-04         1.128e-05              5.3%            1899       2.670085e-06
variance ratio, measured    = 85.8
variance ratio, theoretical = 83.0

What the code does. somma_sui_percorsi is a single function used for two different purposes, and it is the part to read carefully. With the true weights (a, 1−a) it sums the probabilities of the paths touching state 4, that is it returns the exact P(A). With the weights (a²/a′, (1−a)²/(1−a′)) it instead sums p²/q, which is by definition E_q[w²·1_A], the second moment of the importance estimator. The coda factor counts the steps remaining after absorption: the simulator does not stop when the event is settled, and the weight keeps multiplying to the end of the horizon. With the true weights coda = a + (1−a) = 1 and the factor disappears; with the squared weights it is 0.02 + 1.62 = 1.64 and it does not. episodi simulates all episodes in parallel with array operations, accumulating the logarithm of the weight so as not to lose precision.

What it tells us. With 20 000 episodes the naive estimator finds four saturations in total. Its estimate (2.0·10⁻⁴) is not even that far from the truth, but its declared relative error is 50% — and that 50% is itself computed on four events, so it is a measurement that cannot be trusted. It is the worst condition an estimator can be in: it is so imprecise that it cannot measure how imprecise it is. The importance estimator, with the same number of episodes, has 1899 useful ones and reaches 5.3%. To match it, the naive one needs two million episodes: a hundred times the work for the same result.

The variance ratio, and why we report it twice. Measured 85.8, theoretical 83.0: a 3.4% gap, within the uncertainty with which a variance can be estimated on 437 useful events. The theoretical value is not an assertion taken from a book, the function above computes it. It is worth saying that a first version of this derivation gave 137.7 — wrong, because it summed the weights only up to the instant of absorption instead of over the whole horizon. The mistake was invisible in the formula and showed up only as disagreement with the measurement: exactly the check you lose if you report a single number.

When Monte Carlo is the wrong choice

No technique in this lesson is universally better, and this is the section courses leave out.

The rule to take away. Monte Carlo is the method you use when the others do not apply, and almost always they do not. But «almost always» is not «always», and the cost of not noticing is a factor of ten million.

Typical mistakes (symptom, cause, remedy)

1. Seed not fixed. Symptom: two runs give different numbers and nobody knows which to report. Remedy: rng = np.random.default_rng(42), the generator passed as an argument, the seed written next to the result.
2. Using NumPy's global state. Symptom: the result changes when a dependency is updated or a line is moved. Cause: np.random.seed() and np.random.rand() share a state with the whole process. Remedy: default_rng, and SeedSequence.spawn for parallel streams.
3. Reporting the standard deviation of the data instead of that of the estimator. Symptom: an error margin that does not shrink as samples increase. Cause: σ̂ was printed instead of σ̂/√n. Remedy: if the number does not fall with n, it is not a standard error.
4. Trusting a confidence interval on a biased estimator. Symptom: the interval covers the true value while the repetitions are few and stops covering it as they grow — the 1600-repetition row of the queue experiment. Cause: the interval measures the error around the estimator's own mean, not the distance from the truth. Remedy: discard a warm-up, do not truncate the long observations, validate on a case with a known answer.
5. Proposal with lighter tails than the target. Symptom: the estimate seems to converge and then jumps; the sample variance grows with n. Cause: ∫ f²p²/q diverges, the estimator has infinite variance while remaining unbiased. Remedy: print the maximum weight and the ESS; when in doubt, heavier tails.
6. Summing importance weights only up to the stopping time. Symptom: theory predicts a variance reduction systematically larger than the one you measure, and the gap does not close as samples grow. Cause: the likelihood ratio keeps accumulating after the event has been settled, and on a single move the sum of the p²/q is not 1. Remedy: sum the weights over the whole horizon actually simulated.
7. Unnormalised weights, or normalised when they should not be. Symptom: the estimate is wrong by a constant factor, the same at every run. Cause: p or q known only up to the constant. Remedy: the raw estimator (1/n)Σwᵢf(xᵢ) is unbiased but needs exact densities; the self-normalised one Σwᵢf(xᵢ)/Σwᵢ tolerates unknown constants but is biased to order 1/n.
8. Using rejection in high dimension. Symptom: the program runs for hours and produces very few samples, with no errors. Cause: 1/M collapses exponentially — 0.25% at ten dimensions, 2.5·10⁻⁸ at twenty. Remedy: below 1% acceptance you change method, not machine.
9. Choosing n «because it seemed enough». Symptom: nobody can say what precision the result has. Remedy: decide the precision e first, estimate σ on a pilot run, set n ≥ (σ/e)². It is the same formula read backwards, and it turns a simulation into a measurement.

The nine mistakes fall into two families. Five concern the declaration of uncertainty (seed, global state, standard error confused with the standard deviation, interval on a biased estimator, n chosen by eye) and are avoided with a check on a case with a known answer. Four concern proposals and weights (tails that are too light, stopping time, normalisation, rejection in high dimension) and are recognised by the indicators the lesson always has you print: maximum weight, ESS, acceptance rate.

Best practice

• Declare the seed and the versions. Without a seed a result is not reproducible; without a version it is not reproducible a year later.
• Report the standard error next to the estimate. A number on its own is not a measurement; if the standard error cannot be computed, use the bootstrap.
• Size n from the error you want, with n ≥ (σ/e)² from a pilot run, instead of picking a power of ten.
• Validate on a case with a known answer before using one without. With no reference you cannot tell a bug from a result.
• Watch the weights: maximum weight, ESS and acceptance rate must always be printed, not only when something looks off.
• Try variance reduction before increasing n. Here a three-line control variate returned 61 times: the same precision by brute force would cost sixty-one times the samples.

The six rules are ordered from the cheapest to the most demanding: seed, standard error and sizing of n cost one line, validation on a known case costs half an hour, variance reduction requires thinking about the estimator. If you are starting from zero, apply them in the order in which they are written.

Minimal glossary

• Unbiased estimator: its expected value is exactly the quantity sought. It guarantees nothing about a single run.
• Standard error: standard deviation of the estimator, σ/√n. It is not σ, which is the spread of the data.
• Inverse transform sampling: X = F⁻¹(U) with U uniform. Exact and with no waste, requires F⁻¹ in closed form.
• Rejection sampling: propose from q, accept with probability p/(Mq). Exact, at a cost of M attempts per sample.
• Proposal: the distribution q you can draw from, used in place of the target p.
• Importance weight: w = p/q, corrects having sampled from the wrong distribution.
• Effective sample size: ESS = (Σwᵢ)²/Σwᵢ², how many equal-weight samples yours are really worth.
• Antithetic variates: pairs (U, 1−U), negative correlation and lower variance on monotone functions.
• Control variate: a correlated function of known mean, subtracted with c* = Cov(f,g)/Var(g); reduces the variance by the factor 1 − ρ².
• Reparametrisation trick: rewrite x ∼ p_θ as x = g_θ(ε) with ε from a fixed distribution, in order to differentiate inside the expectation.

The glossary has a hidden order: estimator, error, three ways of sampling, variance reduction, reparametrisation. The entries that get confused most often are the first two: an unbiased estimator is not a precise estimator, and the standard error is not the standard deviation of the data.

Mini-project

One script, one question: what is the 99th-percentile response time of a service made of three calls in series?

1. Simulate three independent latencies, each generated by inverse transform from an exponential with a different λ (0.2, 0.5 and 1.0 ms⁻¹), and sum them. Estimate the mean and the 99th percentile, declaring the standard error of both (the percentile needs the bootstrap: there is no closed form);
2. decide n starting from the precision: you want the 99th percentile to better than 0.1 ms. Estimate σ with a pilot run of 10 000, compute the required n, check afterwards that you hit it;
3. add a timeout T beyond which the call fails, and estimate the failure probability for large T, first with naive Monte Carlo and then by importance (hint: propose the exponentials with a smaller λ and weight by the ratio of the densities). Report the two relative errors and the useful episodes of each.

The three requests build a small, complete measurement project: estimate with error, sizing of the samples starting from the precision, rare event with importance sampling. If you cannot state the standard error of the 99th percentile, you have found the point where the bootstrap is needed.

Deliverable: the script must print, for each point, estimate, standard error and number of samples, plus a line saying whether the precision required at point 2 was reached.

Challenge

• Easy. Pen and paper: you have an estimator with standard error 0.08 on n = 5000 samples. How many are needed to bring it to 0.04, and why? Solution: the error falls as 1/√n, so halving it requires quadrupling n: you need 20 000 samples. In general n_new = n · (e_old/e_new)², and the exponent 2 is the one thing to remember from the whole lesson.
• Medium. Pen and paper: derive the inverse transform for the exponential of parameter λ starting from F(x) = 1 − e^{−λx}, and for a discrete variable with three outcomes of probabilities 0.5, 0.3 and 0.2. For the discrete one, write out the three intervals of u explicitly and say what the formula returns for u = 0.5 exactly, justifying the choice with the definition of generalised inverse.
• Hard. Code: choose a density with a known maximum (a different Beta, or a truncated mixture of Gaussians) and implement rejection with a uniform proposal. Measure the rate over 10⁵ attempts against 1/M; then redo it all with M quadrupled and verify that the samples have the same distribution (three quantiles) but cost four times as much.
• Expert. Code: study how the gap between score function and reparametrisation depends on μ. The theoretical variances are μ⁴ + 14μ² + 15 and a constant 4σ²: plot the ratio for μ from 0 to 5, verify it numerically, then add the optimal baseline and measure how much of the gap it recovers as μ grows.
• Enterprise. Take a quantity from your own work that today is reported without uncertainty — a mean response time, a conversion rate, the expected cost of a batch — and estimate it by simulation. Declare the standard error, justify the number of samples by the precision required rather than by the time available, and write one line on what bias there might be in the way the data was collected.

The five levels climb from the computation by hand to the estimate on your own data: the first two consolidate the standard error and the inverse transform, the third and the fourth require code and comparisons with the theoretical value, the last one brings the declaration of uncertainty onto a quantity that you report today without an error. Do first the one where you get stuck.

The code exercises of this lesson are in a single file, esercizi.py (about 15 KB). Three code exercises on domains new with respect to the article: inverse transform for a heavy-tailed Pareto, where the standard error exists as a number but means nothing; the coupon collector, with the exact value n·H_n to be recovered; a control variate with the optimal coefficient computed by hand. A skeleton with the parts to complete and, at the end, the commented solutions. It requires Python 3.12+ and NumPy 2.3. Run it with python3 esercizi.py: it works even with the parts to complete left empty and tells you which ones are missing.

Recap

Mental map. Î_n = (1/n)Σ f(xᵢ) estimates E_p[f], is unbiased by linearity, and has standard error σ/√n — independent of dimension, which is the whole point. Inverse transform: F⁻¹(U), exact but you need the inverse. Rejection: p ≤ Mq, acceptance 1/M, dies beyond five dimensions. Importance: E_p[f] = E_q[f·p/q], variance ∫f²p²/q − I², blows up if q has light tails. Antithetic variables and control variates lower σ instead of raising n. Reparametrisation: x = μ + σε moves the randomness out of the gradient path.

Cheat sheet.

3 concepts: the unbiased estimator with error σ/√n independent of dimension; importance sampling as an exact rewriting that moves samples to where the information is; reparametrisation as moving the randomness out of the gradient path. 3 mistakes: confusing the spread of the data with that of the estimator; trusting a confidence interval on a biased estimator; using a proposal with lighter tails than the target. 3 best practices: declaring seed and versions; sizing n from the error you want; validating on a case with an exact answer before using one without. 1 exercise: the mini-project of this lesson.

Next lesson: Markov chains, that is the processes that forget the past — and the distribution that in the end stops changing.

Questions to take further

• If rejection dies in high dimension and importance requires a good proposal, how do you sample from a complicated distribution in a hundred dimensions? (Markov chain Monte Carlo: lesson 15 of this module for the theory, and Metropolis-Hastings as the application.)
• Where exactly does the stationary distribution used here as an exact reference come from? (It is the left eigenvector with eigenvalue 1 of the transition matrix: the lesson on eigenvalues and SVD provides the tool, the next one applies it.)
• If the estimated mean has a standard error, how large must the difference between two estimates be before you can call them different? (That is the hypothesis-testing problem: the lesson on statistics frames the question, here we have seen only the precision side of it.)

The three questions connect the lesson to the others: the first leads to Markov chain Monte Carlo as the remedy for rejection in high dimension, the second to the stationary distribution as the left eigenvector of the transition matrix, the third to hypothesis testing, that is the comparison between two estimates that both have a standard error. They do not need to be solved right away.

The same series also includes: 28 - Bayesian Statistics: Priors, Posteriors, Conjugates and Why Regularisation Is a Belief (the previous lesson), 30 - Markov Chains: The Property That Forgets the Past, and the Distribution That Stops Changing (the next lesson).