The Gold Standard of Scientific Computing | Monte Carlo Standard Error, Sample Size and Correlated Draws

Suppose a model produces a binary outcome Y, with Y=1 when an event occurs and Y=0 otherwise. Its event probability is p=E[Y]. Run the simulator n times under a defined random mechanism and estimate p using p_hat=(Y₁+…+Yₙ)/n.

If the probability is unknown, p_hat is an estimate, not a newly discovered deterministic property that must be exactly equal to p. Its value changes between independently repeated batches of simulations. The purpose of Monte Carlo error analysis is to quantify that variation under the sampling mechanism.

For a teaching example, imagine the model’s true probability is p=0.30, known to us only so calculations can be checked. A sample of 100 independent trials might produce 30 successes. Then p_hat=0.30, perfectly matching p by coincidence; this does not mean the estimator has zero sampling uncertainty.

A real scientific simulation may estimate a mean, integral, risk or other statistic rather than a binary probability. The Bernoulli example is deliberately simple because its variance and Monte Carlo standard error can be derived exactly under independence.

Contents · Next · Subject root

A single Bernoulli outcome is zero or one. Its variance is Var(Y)=p(1−p). At p=0.30 this is 0.21, so the standard deviation of one outcome is √0.21≈0.4583, or 45.83 percentage points on a zero-to-one scale.

That is not the uncertainty of the estimated probability after many runs. For n independent identically distributed Bernoulli draws, averaging divides the variance by n: Var(p_hat)=p(1−p)/n. The corresponding Monte Carlo standard error is SE(p_hat)=√[p(1−p)/n].

At n=100 the standard error is √(0.21/100)≈0.04583, or about 4.58 percentage points. It is smaller than the variability of individual outcomes because the average combines many independent observations.

These two standard deviations answer different questions: how variable is one model outcome, and how variable is the sample-based estimate across repeated Monte Carlo experiments? A report that calls the raw event standard deviation its estimator precision confuses the two.

Contents · Previous · Next · Subject root

Assume exactly 30 successes and 70 failures occurred among 100 independent simulations. Then p_hat=30/100=0.30. An estimated standard error using the plug-in probability is √[0.30×0.70/100]=0.04583.

Under the hypothetical known-true p=0.30 model, this coincides numerically with the theoretical standard error, but that equality depends on our constructed batch. Another batch with 25 successes would have p_hat=0.25 and a plug-in standard error √(0.25×0.75/100)≈0.04330.

The change in the plug-in standard error does not mean the simulation became more informative simply because fewer events happened; the underlying Bernoulli variance depends on the actual p, which is generally unknown. When probabilities are very close to zero or one, normal approximations and plug-in estimates need further care.

A careful results table records number of runs, number of events, estimated probability, standard-error method and any dependence or variance-reduction techniques used. Printing ‘30.000%’ alone hides the experiment’s limited precision.

Contents · Previous · Next · Subject root

For fixed independent Bernoulli probability p, SE(p_hat) is proportional to 1/√n. Thus multiplying n by four cuts Monte Carlo standard error in half; multiplying n by a hundred cuts it by a factor of ten.

With p=0.30, n=100 gives SE≈0.04583, n=400 gives ≈0.02291, n=1000 gives ≈0.01449 and n=10000 gives ≈0.00458. These correspond to roughly 4.58, 2.29, 1.45 and 0.46 percentage points.

Why not claim four hundred runs are ‘four times more accurate’ than one hundred? Accuracy can refer to bias, mean squared error or other metrics. The standard error falls by two, while the variance falls by four. The metric must be named.

The cost can be significant: cutting the random error by a factor of five ordinarily requires around 25 times as many independent runs if no other method changes. Better algorithms and variance reduction may sometimes improve efficiency, but sample count alone follows this square-root law.

Contents · Previous · Next · Subject root

For an event probability not too near zero or one and a sufficiently large independent sample, the central limit theorem motivates an approximate interval p_hat±1.96 SE_hat for a nominal 95% confidence procedure. For p_hat=0.30 and n=100 this gives approximately 0.30±0.0898, or about 21.0% to 39.0%.

This is an approximate statistical interval, not a guarantee that the true p is in a particular realized interval with 95% probability under a frequentist interpretation. The method is designed to cover the fixed unknown parameter in about 95% of repeated eligible experiments under its assumptions.

At small n or extreme p_hat, the simple symmetric normal interval can perform poorly and may even extend below zero or above one. Better-calibrated methods, such as Wilson intervals or exact binomial intervals, may be preferable for those cases.

Do not confuse this Monte Carlo sampling interval with physical model uncertainty. A model could estimate a probability very precisely yet represent the real world poorly because of incorrect inputs, structural assumptions or omitted processes.

Contents · Previous · Next · Subject root

Suppose you run ten independent batches of 100 simulations. Each batch produces a probability estimate, and the empirical variation among those ten estimates offers a rough check of Monte Carlo uncertainty. The resulting batch estimates are not ten new model parameters; they are repeated measurements of the same underlying probability.

Pooling all independent outcomes into one batch of 1000 draws generally provides a more precise estimate than quoting one isolated batch of 100. Under p=0.30 the ideal standard error for 1000 draws is about 0.01449.

However, batch-level variation is itself noisy with only ten batches. If the code accidentally reuses the same random seed and regenerates identical output each time, the empirical batch variance may be zero even though the sampling uncertainty is not.

Random seeds are part of reproducibility, not a substitute for independence. Recording a seed helps reconstruct a run; choosing independent streams appropriately is necessary if the study relies on independence across replications.

Contents · Previous · Next · Subject root

For draws Y₁,…,Yₙ that may be dependent, Var(mean)=(1/n²)ΣᵢΣⱼCov(Yᵢ,Yⱼ). Under independence all cross-covariances vanish, leaving Var(Y)/n. When nearby draws are positively correlated, those extra covariance terms can enlarge estimator variance.

For an approximately stationary chain with common marginal variance σ² and autocorrelation ρ_k at lag k, a long-run approximation is Var(mean)≈(σ²/n)[1+2Σ from k≥1 ρ_k], provided the sum behaves appropriately.

The bracketed factor describes how serial dependence changes the information carried by n draws. It leads to an effective sample size n_eff≈n/[1+2Σρ_k]. Positive autocorrelation usually reduces effective sample size; some anticorrelated methods can increase it under suitable conditions.

This is why ‘one million iterations’ is not a complete uncertainty statement for a Markov chain Monte Carlo calculation or a correlated time-series simulation. The dependence structure and mixing behaviour are essential.

Contents · Previous · Next · Subject root

Construct an idealised stationary sequence whose lag-k autocorrelation follows ρ_k=0.8^k. The sum of positive lags is 0.8/(1−0.8)=4, so the long-run variance inflation factor is 1+2×4=9.

Therefore a long-run effective sample-size approximation is n_eff≈n/9. For n=100, this gives about 11.1 independent-equivalent draws, rather than 100. At p=0.30 the corresponding standard-error approximation is √(0.21/11.1)≈0.1375, much larger than 0.04583.

The n/9 calculation is asymptotic, because a finite record has an autocorrelation-weighted sum with factors depending on lag. It should be presented as an illustrative long-run approximation, not an exact finite-n confidence interval.

The point is conceptual but consequential: collecting many neighbouring states from a sticky simulation can yield less information than fewer independent draws. Longer runs, better transitions or suitable variance reduction may be needed.

Contents · Previous · Next · Subject root

Effective sample size is estimated from a finite correlated chain; it is not usually known as exactly as n. Empirical autocorrelations at large lags are noisy, and different estimators can truncate or regularise their sums differently.

For multiple Markov chains, separate concerns include warm-up, convergence diagnostics, mixing, multimodality and whether the chain reached the distribution being targeted. A large reported effective sample size cannot rescue a chain trapped in the wrong region of parameter space.

Stan’s reference manual describes effective sample size in terms of autocorrelation and Monte Carlo standard error. That framework is useful, but a reader should still ask which quantity’s ESS was calculated; different parameters or functions of one simulated state can have different correlations.

One global ESS for all results is often an oversimplification. Report uncertainty for the actual estimated statistic and explain how dependence was assessed.

Contents · Previous · Next · Subject root

Monte Carlo standard error concerns variation due to finite random sampling under a specified simulation procedure. It does not automatically account for programming mistakes, discretisation error, measurement bias, mis-specified model parameters or systematic approximation in the simulation algorithm.

For example, doubling a simulation run from 10,000 to 20,000 draws reduces ideal random error by about √2, but if the simulator uses the wrong physical rate constant, its estimate may converge very precisely to the wrong value.

Likewise, a pseudorandom generator with flawed implementation or an algorithm that never explores part of its intended distribution can invalidate the independent-draw formulas. More iterations do not correct a missing model state.

A trustworthy computational report separates sampling precision, numerical approximation and model validity. A small standard error is one component of evidence, not the whole meaning of accuracy.

Contents · Previous · Next · Subject root

Exercise A: A binary event has p=0.30 and n=400 independent simulations. Its theoretical standard error is √(0.21/400)≈0.02291, or 2.29 percentage points. It is half the standard error from n=100.

Exercise B: To reduce the ideal standard error from about 0.0458 to about 0.00458 without changing p or estimator, multiply the independent draw count by 100. One hundred runs become ten thousand.

Exercise C: A stationary chain has ρ_k=0.8^k. What is the long-run variance inflation factor? It is nine, yielding effective size about n/9. Explain why counting adjacent states as independent exaggerates precision.

Exercise D: A colleague reports 30% with zero uncertainty because three batches using the same seed all returned identical results. The batches are not evidence of independent replication. The seed and random-stream design must be examined.

Contents · Previous · Next · Subject root

Identify the target statistic, sampling algorithm and assumptions about independent or dependent draws. Record n, number of batches, random generator and seed handling, and whether any observations were excluded or selected after seeing results.

For independent binary draws, report p_hat and a suitable standard error or confidence interval; for correlated chains, report estimated effective sample size and an appropriate Monte Carlo standard error with convergence checks. Avoid confusing the outcome standard deviation with uncertainty in the estimated mean.

Check results against a case with known analytic probability before applying the method to an unfamiliar system. Then vary simulation count and, where practical, compare independent repeated runs. A correct computational experiment should show error shrinking at a plausible rate rather than only producing more digits.

Continue from The Gold Standard of Scientific Computing to compensated summation and Euler stability. These distinguish randomness, floating-point effects and discretisation errors.

Contents · Previous · Subject root