Monte Carlo Accuracy Is About Variance, Not Just Samples

Monte Carlo methods are often introduced with a reassuringly simple recipe: generate random draws, evaluate a function, average the results, and increase the sample size until the answer stabilises. The mathematical justification is equally familiar. If independent observations have finite variance, the central limit theorem gives an error that decreases at the rate $N^{-1/2}$. This rate is dimension-independent in a narrow but important sense, and it explains why Monte Carlo remains useful when deterministic quadrature becomes difficult.

The same presentation can encourage a misleading computational habit. If the estimate is noisy, run more simulations. If the desired standard error is ten times smaller, run one hundred times as many. That conclusion follows only after the estimator itself has been fixed. Monte Carlo methods are not merely about how many random numbers are generated. They are about constructing a random estimator of a deterministic quantity, and two unbiased estimators of exactly the same target can have variances that differ by factors of ten, one thousand, or much more.

Variance is therefore a computational resource. Reducing it changes the amount of simulation required for a given precision. Control variates, antithetic sampling, stratification, conditional expectation and importance sampling all exploit known structure to change variance without changing the target quantity. In difficult rare-event problems, estimator design can be the difference between a calculation that finishes in seconds and one that is effectively impossible.

This article concerns classical Monte Carlo integration rather than Markov chain Monte Carlo. The blog already has separate material on MCMC and importance sampling. The purpose here is to develop the common mathematical principle behind variance reduction and to quantify what the gains actually mean.

The square-root law is a statement about one estimator

Suppose the target is

If

are independent draws from $p$, the ordinary Monte Carlo estimator is

Provided

the estimator is unbiased,

with variance

The central limit theorem gives

A natural Monte Carlo standard error is therefore

where $s_f^2$ is the sample variance of the simulated values $f(X_i)$.

The square-root rate is both powerful and unforgiving. If the standard error must be divided by two, the sample size must be multiplied by four. Dividing error by ten requires roughly one hundred times as many draws. Dividing it by one hundred requires roughly ten thousand times as many draws. Brute-force Monte Carlo becomes expensive quickly when each function evaluation involves a PDE solve, a stochastic simulator, a financial portfolio revaluation, a reliability model or another expensive forward computation.

The rate also hides the constant

Two estimators can both converge at rate

while one has a variance one hundred times smaller. At equal cost, its standard error is ten times smaller. Equivalently, the high-variance estimator needs one hundred times as many effective evaluations to match it.

For this reason, the useful comparison is not simply asymptotic rate. It is variance at a fixed computational budget, or more generally a cost-normalised measure such as

A variance-reduction method that halves variance while doubling evaluation cost has achieved nothing at equal budget. A method that reduces variance by a factor of thirty while using the same number of function evaluations has changed the computational problem substantially.

Control variates use a quantity whose expectation is already known

The cleanest variance-reduction argument begins with correlation. Suppose $Y=f(X)$ has unknown mean $\mu$, while another random variable $Z$, generated from the same draw, has known expectation

For any constant $c$,

remains unbiased because

Its single-draw variance is

Minimising with respect to $c$ gives

At this value,

where $\rho_{YZ}$ is the correlation between $Y$ and $Z$. A control variate therefore converts correlation with a known quantity directly into variance reduction.

An exact Gaussian example makes the gain visible. Let

and suppose the target is

The naive Monte Carlo variable is

Using Gaussian moments,

so

A natural control variate is

because

is known exactly. Its variance is

and

The optimal coefficient is therefore

The controlled variable is

and its variance is

The estimator has exactly the same expectation as the naive estimator and one quarter of its variance. With $N$ draws,

while

The naive method therefore needs four times as many samples to obtain the same variance.

This example is modest compared with many practical applications because the correlation is not extremely close to one. In simulation models, useful control variates often come from simplified analytic approximations, conservation laws, lower-fidelity simulators or quantities whose expectation is available by symmetry. When a control tracks most of the variation in the expensive output, the gain can be dramatic.

The central principle is that randomness in a Monte Carlo estimator is not automatically useful. If part of that randomness is shared with a quantity whose mean is already known, it can be subtracted rather than averaged away slowly through larger $N$.

Deliberate dependence can reduce variance

Independent sampling is convenient, not sacred. Negative dependence can make an average more stable.

Let

and consider

One ordinary estimator evaluates $e^U$ at independent uniforms. Antithetic sampling uses each draw together with its reflection,

The pair estimator is

Because both $U$ and $1-U$ are uniformly distributed,

The two function evaluations are negatively correlated. When $U$ is large, $1-U$ is small, and vice versa.

For one ordinary evaluation,

For one antithetic pair,

A fair comparison should use the same number of function evaluations. With $2N$ independent evaluations, ordinary Monte Carlo has variance

With $N$ antithetic pairs, also requiring $2N$ evaluations, the variance is

The variance ratio is therefore approximately

At equal evaluation cost, antithetic coupling is more than thirty times as efficient in this example.

There is no universal guarantee that antithetic sampling helps. The construction works when the two evaluations tend to move in opposite directions around the target. For monotone functions of a scalar uniform variable, the reflection $U\mapsto1-U$ often induces useful negative correlation. For irregular or non-monotone functions, the covariance may be weak or even positive. The right diagnostic is the covariance of the paired evaluations, not the fact that the method is called antithetic.

Stratified sampling exploits a related idea by forcing the sample to cover parts of the domain deliberately rather than trusting random allocation. If the integration domain is partitioned into strata $S_h$ with probabilities $p_h$, then

Estimating each conditional expectation separately prevents random overrepresentation of some strata and underrepresentation of others. With suitable allocation, variance can fall substantially when outcomes are relatively homogeneous within strata.

Conditional Monte Carlo goes further. If $Y$ is a noisy estimator and $Z$ is a useful conditioning variable, replace $Y$ by

whenever that conditional expectation can be computed. The law of total variance gives

Therefore,

This is the variance-reduction version of Rao-Blackwellisation. Randomness that can be integrated out analytically should not be simulated merely because a simulator can generate it.

Rare events expose the limits of brute-force simulation

Variance reduction becomes essential when the target itself is rare. Consider

and the probability

The exact value is approximately

A naive Monte Carlo estimator is

This is simply a Bernoulli average. Its variance is

For rare events, absolute standard error is less informative than relative error. The relative standard error is approximately

To obtain relative standard error near 10%, we need roughly

so

That is about one hundred billion independent standard-normal draws for an estimator whose relative uncertainty is still around ten percent.

The problem is not that normal random numbers are expensive. The problem is that almost every simulated observation contributes exactly zero. At

the expected number of six-sigma exceedances is only about

The probability of observing no exceedance at all is approximately

A simulation containing one hundred million draws will therefore return the estimate zero about ninety percent of the time.

Importance sampling changes the distribution from which simulations are drawn. Let the target density be the standard-normal density $\phi(x)$, but simulate instead from

the density of

The probability can be written exactly as

The likelihood ratio is

The importance-sampling estimator is therefore

with

The estimator remains unbiased, but the event now occurs with probability one half under the proposal rather than probability $10^{-9}$.

Its second moment can be calculated analytically:

Since

the relative variance of one importance-sampling draw is

The relative standard error of the sample mean is therefore

For 10% relative error, this requires roughly

draws.

The contrast is not subtle.

Method	Approximate draws for 10% relative standard error
Naive Monte Carlo	$1.01\times10^{11}$
Shifted importance sampling $N(6,1)$	$687$

Both estimators target exactly the same probability. Both are unbiased. Both ultimately obey a square-root law. The difference lies in the variance constant.

This is why statements such as "Monte Carlo converges at $N^{-1/2}$" are mathematically correct and computationally incomplete. The exponent of $N$ does not tell us whether the constant multiplying it is $10^{-3}$, $1$, or $10^9$.

The zero-variance importance distribution would sample directly from the target distribution conditional on the rare event,

If this proposal were available, every weighted sample would equal the desired probability and the estimator would have zero variance. Of course, constructing the exact conditional law usually requires knowledge as difficult as the original problem. The theoretical optimum is still useful because it reveals what a good practical proposal should imitate: it should concentrate probability where the integrand contributes most.

Importance sampling can also be much worse than naive Monte Carlo

Changing measure does not automatically improve an estimator. For

importance sampling from $q$ uses

Its variance is finite only if

This condition explains why proposal tails matter. If $q(x)$ becomes too small in regions where $f(x)p(x)$ remains important, the importance weights

can become enormous. A few observations then dominate the estimate. The method can have much higher variance than naive Monte Carlo and, in extreme cases, infinite variance even though the target expectation itself is finite.

Support is an even more basic requirement. If

on a region where

then the proposal never visits part of the integral. No reweighting can recover a region that is assigned zero sampling probability.

This is why stable importance sampling should be assessed through its weights as well as its final estimate. A common descriptive diagnostic is the effective sample size,

or, for normalised weights $\tilde w_i$,

This quantity is only a heuristic summary and should not be confused with a universal variance identity. It is nevertheless useful for exposing weight concentration. Ten thousand simulated points can carry the information of only a few observations when almost all total weight is assigned to a tiny subset.

Self-normalised importance sampling,

is useful when the target density is known only up to a normalising constant. Unlike the basic importance-sampling estimator, it is generally biased at finite $N$, although it can be consistent under appropriate conditions. The distinction illustrates a broader principle: Monte Carlo design is not limited to choosing between unbiased estimators. A small controlled bias can sometimes be worthwhile if it reduces variance or makes the computation possible, but the trade-off should be stated explicitly.

Variance reduction should be evaluated at equal cost

The previous examples counted function evaluations because that is the dominant cost in many simulations. Real Monte Carlo algorithms have more complicated cost structures. A control variate may require an additional model evaluation. A stratified design may require preprocessing. Importance sampling may need optimisation to build the proposal. Conditional Monte Carlo may replace simulation with an expensive numerical integral.

The relevant comparison is therefore not variance per draw but variance per unit of computation. If method $A$ costs $c_A$ per independent replication and has single-replication variance $v_A$, while method $B$ has cost $c_B$ and variance $v_B$, then at a total budget $C$,

and similarly for $B$. The product

is a natural first-order measure of inefficiency.

This matters in multi-fidelity simulation. Suppose a high-fidelity model is expensive and a lower-fidelity approximation is cheap but highly correlated with it. The cheap model can be used as a control variate. A large number of low-fidelity evaluations can then reduce the variance of a much smaller number of high-fidelity simulations. The optimal allocation depends on both correlations and relative costs.

The same logic appears in nested Monte Carlo, where an outer simulation calls an inner simulation. Increasing inner accuracy uniformly can be spectacularly wasteful if many outer states contribute little to the final quantity. Adaptive allocation can devote computation where conditional variance is large or where the decision boundary is uncertain.

Monte Carlo computation should therefore be budgeted around variance contribution, not around a fixed number of samples at every stage.

Standard errors are part of the result

A simulation estimate without a Monte Carlo error assessment is incomplete. For independent finite-variance draws, the estimated MCSE

is straightforward. A simulation can stop when the MCSE falls below a scientifically meaningful absolute or relative tolerance rather than after an arbitrary round number of iterations.

The stopping criterion should be connected to the scale of the decision. If the estimate is approximately $1000$, an MCSE of $0.01$ may be pointless if model uncertainty is $50$. Conversely, a simulation standard error of $10$ may be unacceptable if the difference between two policies is expected to be $5$.

Monte Carlo uncertainty is not the same as uncertainty in the underlying scientific model. If a posterior mean is estimated by simulation, the posterior uncertainty describes uncertainty about the model parameter conditional on the model and data, while the MCSE describes numerical uncertainty from approximating that posterior expectation with finite simulation. These two layers should not be conflated.

Variance estimation also becomes more difficult when draws are dependent. In Markov chain Monte Carlo,

are intentionally correlated. For a stationary chain, the asymptotic variance of the sample mean can often be written as

where

is an integrated autocorrelation time. Positive serial correlation reduces the effective information in $N$ iterations. This motivates the familiar approximation

The dependent case reinforces the main theme rather than changing it. Raw iteration count is not a measure of accuracy. Variance is.

The square-root law has assumptions

The classical central limit theorem requires finite second moments and appropriate independence or weak-dependence conditions. Monte Carlo estimators can violate these assumptions.

If

the usual

standard error need not apply. Heavy-tailed simulation outputs can produce estimates dominated by rare enormous contributions, with sample variance itself unstable. Running longer may create the appearance of convergence for a while and then change the estimate dramatically when a previously unseen tail event occurs.

This is especially relevant in importance sampling because poor weight distributions can create infinite second moments even when the original integrand has finite variance. Weight diagnostics, tail analysis and proposal sensitivity are therefore part of numerical reliability, not optional extras.

High dimension creates a different limitation. Plain Monte Carlo's nominal convergence rate does not contain dimension explicitly, which is often contrasted favourably with deterministic grids. The variance constant can still deteriorate badly with dimension, and importance proposals can become mismatched in exponentially many directions. A dimension-free exponent does not imply dimension-free computational difficulty.

Quasi-Monte Carlo methods attack integration error differently by replacing independent random points with low-discrepancy sequences. Under smoothness and effective-dimension conditions, they can outperform the classical square-root rate. Randomised quasi-Monte Carlo restores a form of error assessment while retaining much of the low-discrepancy advantage. This is a distinct subject deserving its own treatment; it is enough here to note that $N^{-1/2}$ is not a universal computational lower bound for numerical integration.

Monte Carlo is an estimator-design problem

The examples lead to a consistent conclusion. Ordinary Monte Carlo estimated

with single-draw variance $96$. A control variate reduced that variance to $24$ without changing the expectation. Antithetic sampling estimated

more than thirty times as efficiently as independent sampling at equal evaluation cost. For the six-sigma Gaussian tail, naive simulation required about $10^{11}$ draws for 10% relative precision, while a simple shifted importance sampler required fewer than $10^3$.

None of these improvements changed the quantity being estimated. They changed the random variable used to estimate it.

That is the most useful way to think about Monte Carlo methods. The problem is not merely to draw random numbers from the most obvious distribution and average them. The problem is to construct an estimator whose expectation is correct, whose variance is controlled, whose computational cost is acceptable, and whose error can be measured.

Sometimes the best method is indeed ordinary sampling. It is simple, parallel, auditable and surprisingly robust. When function evaluations are cheap and the variance is moderate, elaborate variance reduction can add complexity without practical benefit.

When simulation is expensive or the target is rare, the estimator itself becomes the main numerical object. Known expectations can become control variates. Symmetry can create antithetic pairs. Conditioning can remove unnecessary randomness. Stratification can prevent inefficient allocation. A change of measure can move computation into a rare region that naive simulation almost never visits.

The $N^{-1/2}$ law remains true for a large class of these estimators. It simply does not tell us how large $N$ needs to be.

That number is controlled by variance.

References

Bucklew, J. A. (2004). Introduction to Rare Event Simulation. Springer.

Glasserman, P. (2004). Monte Carlo Methods in Financial Engineering. Springer.

Hammersley, J. M., & Handscomb, D. C. (1964). Monte Carlo Methods. Methuen.

Kahn, H., & Marshall, A. W. (1953). Methods of reducing sample size in Monte Carlo computations. Journal of the Operations Research Society of America, 1(5), 263–278.

Liu, J. S. (2001). Monte Carlo Strategies in Scientific Computing. Springer.

Owen, A. B. (2013). Monte Carlo Theory, Methods and Examples.

Robert, C. P., & Casella, G. (2004). Monte Carlo Statistical Methods (2nd ed.). Springer.

Rubinstein, R. Y., & Kroese, D. P. (2016). Simulation and the Monte Carlo Method (3rd ed.). Wiley.