Why Monte Carlo Error Falls as the Square Root of the Sample Count

An estimator that averages $n$ independent draws converges, but it converges slowly and without a derivative. Each extra decimal digit of accuracy costs a hundred times the work, and the cost does not depend on the dimension of the problem or on the smoothness of the integrand. That indifference is the method’s selling point and its limitation: it is the only general technique whose rate survives a hundred-dimensional integral, and it is a poor choice for a smooth one-dimensional integral where a quadrature rule converges far faster.

The variance of a mean, and nothing more

For independent identically distributed draws $X_1, \dots, X_n$ with mean $\mu$ and variance $\sigma^2$, the estimator $\hat\mu = \frac{1}{n}\sum X_i$ has

$$\mathrm{Var}(\hat\mu) = \frac{\sigma^2}{n}, \qquad \text{standard error} = \frac{\sigma}{\sqrt{n}}.$$

Independence is what turns the variance of the sum into $n$ times the variance of one term. Drop it and the cross terms appear: for correlation $\rho$ between successive draws, the variance of the mean becomes

$$\frac{\sigma^2}{n}\left(1 + 2\sum_{k=1}^{n-1}\left(1 - \frac{k}{n}\right)\rho_k\right),$$

which converges to $\sigma^2 \tau / n$ where $\tau = 1 + 2\sum_{k \ge 1} \rho_k$ is the integrated autocorrelation time. For an MCMC sampler $\tau$ is routinely in the tens, so a million draws buy the error of forty thousand independent ones. The square-root law is unchanged; the constant is what the sampler’s memory of its own past multiplies.

The $\sqrt{n}$ is not an artifact of the estimator. For any unbiased estimator, the Cramér-Rao bound ties the variance to the Fisher information, and when the information grows linearly in $n$ the standard deviation can only shrink as $1/\sqrt{n}$. Randomised quadrature, quasi-Monte Carlo with low-discrepancy points, and stratified designs beat the constant rather than the exponent: they reduce the effective variance, and their advantage over plain sampling narrows as the dimension rises.

Reading the rate as a work budget

To halve the standard error, quadruple the sample. To gain one decimal digit, multiply the sample by a hundred. To gain two digits, multiply by ten thousand. The progression is worth keeping in view when a tolerance is specified: a relative error of $10^{-3}$ needs a variance-to-mean-squared ratio of about $10^6/n$, so an estimator with a coefficient of variation near 1 needs roughly a million draws, and one with a coefficient of variation near 10 needs a hundred million.

Dimension enters only through the variance. Squeezing a function into a unit hypercube does not increase $\sigma^2$, though the variance of a badly scaled integrand can be enormous: an indicator function whose probability is $p$ has variance $p(1-p)$, which is at most $1/4$, but a function with a rare large contribution has a variance dominated by the rare event and a coefficient of variation like $1/\sqrt{p}$. This is why the same nominal sample size gives wildly different precision for different integrands, and why a variance estimate is not optional.

There is one structural advantage the rate conceals. Quadrature error for a $d$-dimensional tensor grid with $m$ points per axis is $O(m^{-k/d})$ for a rule of order $k$, so a fixed error target demands a number of points exponential in $d$. Monte Carlo’s $n^{-1/2}$ has no $d$ in the exponent at all. The crossover is around the dimension where the quadrature constant’s dimension penalty outweighs the rate advantage, which for smooth integrands is typically somewhere in the single digits to low tens.

Estimating the error bar from the sample

$\sigma$ is rarely known, so the standard error is estimated by replacing it with the sample standard deviation

$$s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(X_i - \hat\mu)^2, \qquad \widehat{\text{se}} = \frac{s}{\sqrt{n}}.$$

The $n-1$ divisor makes $s^2$ unbiased for $\sigma^2$. Using $n$ instead is a bias of order $1/n$ in the variance, which matters only for small samples but is a common and avoidable mistake. The interval $\hat\mu \pm 1.96,\widehat{\text{se}}$ has nominal $95%$ coverage by the central limit theorem, and the coverage is asymptotic: for a skewed $X$ the actual coverage at $n = 100$ can be well below $95%$ while the nominal figure still prints.

Two refinements are standard. Student’s $t$ with $n-1$ degrees of freedom widens the interval to account for the estimation of $s$, which is a visible correction below a few hundred samples and vanishes above a few thousand. The bootstrap resamples the $X_i$ with replacement, recomputes the mean on each resample, and reports the empirical quantiles of those means; it handles skewed and bounded estimators where the normal approximation is poor, at the cost of a second Monte Carlo layer whose own error must be budgeted.

import math, random

def mean_and_error(xs):
    n = len(xs)
    mean = sum(xs) / n
    var = sum((x - mean) ** 2 for x in xs) / (n - 1)
    se = math.sqrt(var / n)
    return mean, se

def batch_means(xs, b):
    # Divide into b batches; scatter of batch means estimates the variance
    # without assuming independence between adjacent draws.
    n = len(xs)
    size = n // b
    sums = [sum(xs[i * size:(i + 1) * size]) / size for i in range(b)]
    grand = sum(sums) / b
    var = sum((s - grand) ** 2 for s in sums) / (b - 1)
    return grand, math.sqrt(var / b)

Batch means is the practical answer for correlated output. Splitting the run into $b$ contiguous blocks and treating the block averages as independent gives a variance estimate that absorbs the autocorrelation, provided each block is long enough that the correlation has decayed. There is no way to check that condition from the data alone, which is why batch-means estimates are usually reported alongside a plot of the running estimate against block length.

Bias, variance, and the mean squared error budget

The square-root law describes the variance of an unbiased estimator. A sampler that is slightly biased trades a term that vanishes as $1/n$ for a term that never vanishes, and the trade is bad once the sample is large. If $\hat\mu$ has bias $b$ and standard deviation $\sigma/\sqrt{n}$, the mean squared error is

$$\mathrm{MSE} = b^2 + \frac{\sigma^2}{n}.$$

Doubling the sample halves the second term and does nothing to the first, so past $n \approx \sigma^2/b^2$ additional work buys no accuracy at all. Every truncation, every discretization of a state space, and every finite-difference approximation inside the estimator contributes to $b$, and the appropriate target for a large run is to drive $b$ below the standard error at the intended sample size, not to drive it to zero.

This is the reason a coarse model is worth checking at a small sample before a fine one is run at a large sample. At $n = 10^3$ a bias of $0.03\sigma$ contributes $9 \times 10^{-4}\sigma^2$ to the MSE, comparable to the variance term $\sigma^2/1000 = 10^{-3}\sigma^2$; at $n = 10^6$ it dominates by a factor of thirty. The bias that is invisible in the pilot run is the one that decides the answer in the production run.

The central limit theorem that justifies the normal interval also needs a moment condition. The classical statement requires a finite variance, and an estimator built on a heavy-tailed quantity can have an infinite variance even when its mean is finite. The sample mean of a Pareto distribution with shape parameter between 1 and 2 converges to the true mean, but the fluctuations around it are of order $n^{1 - 1/\alpha}$ rather than $n^{-1/2}$, and the sample variance itself diverges. An error bar computed from such a sample is meaningless rather than merely imprecise: it does not shrink, and its reported value depends on whichever extreme observations happened to appear.

The rate does not hold for the estimate of the error

The standard error is itself a random variable, and its own relative error is $1/\sqrt{2(n-1)}$ under normality, from the chi-square distribution of $s^2$. At $n = 10^4$ that is about $0.7%$, comfortable. At $n = 100$ it is about $7%$, which means a reported error bar of $\pm 3$ could plausibly be $\pm 3.2$ or $\pm 2.8$ from the same data. Reporting four significant digits on a standard error computed from a hundred draws is a claim about precision that the sample does not support.

Worse, the estimator can converge to the wrong answer while its own error bar shrinks. A sampler with a truncated tail, a generator with lattice structure aligned to the integrand, or a numerical overflow in an intermediate quantity all produce output with a small sample variance around a biased mean. The error bar measures dispersion, not bias, and no amount of sampling detects a bias that is constant across draws. Checking against a known case, a second independent method, or an analytically solvable limit is the only defence.

The same logic explains why two runs that disagree by two standard errors are not necessarily inconsistent. The difference of two independent means has variance equal to the sum of their variances, so the standard error of the difference is $\sqrt{\text{se}_1^2 + \text{se}_2^2}$, larger than either. Comparing runs at the level of their individual error bars overstates the agreement.

Where the square root is not the whole story

Variance reduction changes the constant in front of $n^{-1/2}$, and the reductions are large enough that they change what is feasible. Antithetic pairing can cancel the leading fluctuation when the integrand is monotone in the uniform, giving a variance ratio that is bounded rather than merely improved. Stratification removes the between-stratum component of the variance exactly, so the residual is the within-stratum variance divided by the total sample. Importance sampling moves the work to where the integrand is large but pays for it with a weight variance that can be worse than the original problem.

Each of these is a change to the estimator, not to the rate. The error still falls as $1/\sqrt{n}$ asymptotically, but the constant can be smaller by orders of magnitude, and the crossover to the asymptotic regime can be far out. A practical reading is that a variance reduction of a factor of a hundred is worth a factor of ten thousand in sample count, which is why a modest amount of analysis before a large run pays for itself.

Stratification shows the arithmetic most clearly. Split the sampling region into $k$ strata with probabilities $w_j$ and conditional variances $\sigma_j^2$. Proportional allocation, with $n w_j$ draws in stratum $j$, gives a variance of $\frac{1}{n}\sum_j w_j \sigma_j^2$, against $\frac{\sigma^2}{n}$ for plain sampling, where $\sigma^2$ is the total variance including the between-stratum part. The saving is exactly the between-stratum variance $\sum_j w_j(\mu_j - \mu)^2$, so stratification is worth most when the stratum means differ and cheap to apply when the strata are easy to define. It never hurts, because the between-stratum term is subtracted and the within-stratum term is what remains. It also constrains the design: the strata must be fixed in advance and each must receive its allotted count, which conflicts with samplers whose draw costs vary, such as rejection loops.

The most common failure of all is treating correlated output as independent. The formula $\sigma/\sqrt{n}$ applied to an MCMC chain, or to a discrete-event replication that has not reached steady state, produces an error bar that is too small by the square root of the integrated autocorrelation time. The detailed balance condition guarantees the chain converges to the right distribution, and says nothing about how many samples that takes; the weight variance in importance sampling plays the same role for a reweighted estimator, where the effective sample size rather than the raw count determines the error.