How Accurate Is a Monte Carlo Estimate?

The estimate is an average

Almost every Monte Carlo method estimates an expected value θ = E[X] by drawing independent copies X₁, …, X_n and averaging them. In Monte Carlo π, X is 4 when a random point lands inside the quarter circle and 0 otherwise; in Monte Carlo integration, X is the scaled function value at a random point; in option pricing, it is a discounted payoff.

The average is unbiased, and its standard deviation — the standard error — is σ/√n, where σ is the standard deviation of a single X. Because σ is unknown, it is replaced by the sample standard deviation s:

By the central limit theorem, θ̂ is approximately normal for large n, so θ̂ ± 1.96 · SE is an approximate 95% confidence interval. The central limit theorem simulation shows how large "large" has to be for skewed inputs.

Computing it in one pass

Storing every sample to compute s afterwards is wasteful, and the textbook shortcut Σx² − n·x̄² loses precision badly when the mean is large compared with the spread. Welford’s update keeps a running mean and a running sum of squared deviations and is numerically stable:

// Running mean and standard error (Welford's algorithm).
function estimate(sample, n, rng) {
  let mean = 0, m2 = 0;
  for (let i = 1; i <= n; i++) {
    const x = sample(rng);
    const d = x - mean;
    mean += d / i;
    m2 += d * (x - mean);
  }
  const se = Math.sqrt(m2 / (n - 1) / n);
  return { mean, se, lo: mean - 1.96 * se, hi: mean + 1.96 * se };
}

// E[e^U] for U ~ Uniform(0, 1); the exact value is e − 1.
estimate((r) => Math.exp(r.uniform()), 10000, createRng(42));

def estimate(sample, n, rng):
    """Running mean and standard error (Welford's algorithm)."""
    mean, m2 = 0.0, 0.0
    for i in range(1, n + 1):
        x = sample(rng)
        d = x - mean
        mean += d / i
        m2 += d * (x - mean)
    se = math.sqrt(m2 / (n - 1) / n)
    return mean, se, (mean - 1.96 * se, mean + 1.96 * se)

estimate(lambda r: math.exp(r.uniform()), 10_000, Rng(42))

With the site’s generator and seed 42, both versions return 1.71267 with a standard error of 0.00495; the interval [1.7030, 1.7224] contains the true value e − 1 ≈ 1.71828. Repeating the experiment with 1,000 different seeds (1,000 samples each), the interval covered the true value in 957 runs — close to the promised 95%.

The 1/√n rule and what it costs

Since SE = σ/√n, halving the error takes four times as many samples, and one more correct decimal digit takes a hundred times as many. Solving for n gives the sample size for a target half-width ε at 95% confidence:

For Monte Carlo π, σ = 4√(p(1 − p)) with p = π/4, about 1.642. Pinning π down to ±0.01 therefore needs roughly (1.96 × 1.642 / 0.01)² ≈ 104,000 points; to ±0.001, about 10 million. Estimate σ from a short pilot run, then size the real one.

The rate does not depend on dimension. A grid-based integration rule that is excellent in one dimension needs exponentially many points as dimensions are added, while the Monte Carlo error stays σ/√n. That is why Monte Carlo dominates high-dimensional problems even though it is slow in low dimensions.

Making σ smaller instead of n larger

Because the error is proportional to σ, reducing the variance of what is averaged is as good as buying more samples. Common techniques:

• Antithetic variates: pair each draw u with 1 − u (or z with −z) so that errors partly cancel. The option pricing model lets you switch this on and compare standard errors.
• Control variates: subtract a correlated quantity whose mean is known exactly.
• Importance sampling: sample more often where the integrand matters and reweight.

When samples are correlated

The formula SE = s/√n assumes independent samples. Two common situations break that:

• Markov chain Monte Carlo. Each Metropolis–Hastings sample is a small move from the previous one. Positive autocorrelation means n correlated samples carry the information of fewer independent ones — the effective sample size, roughly n/(1 + 2Σρ_k) where ρ_k is the lag-k autocorrelation.
• Time averages from one long simulation. Consecutive queue lengths in an M/M/1 queue are strongly correlated, especially near full utilization. Split the run into long batches, average each batch, and compute the standard error from the batch means (the method of batch means), or run independent replications with different seeds.

In both cases the naive standard error is too small, and intervals built from it are overconfident.

Bias is a separate problem

A standard error measures random scatter, not systematic error. A time step that is too large in an ODE model, a warm-up period that was not discarded from a queue simulation, or a flawed random-number generator all shift the estimate in a way no amount of sampling removes. Check models against exact results where they exist — every simulation page on this site reports such a check.