Quantifying Monte Carlo Error in Simulation Results

Why quantify Monte Carlo error?

A trial simulation estimates operating characteristics by repeating a specified simulation procedure. For example, each replicate may generate a virtual trial and apply a prespecified decision rule; the replicate is recorded as a success if the rule is met. The fraction of successful replicates estimates the success probability, but it is itself subject to random variation because only a finite number of replicates were run.

Monte Carlo error describes this finite-replicate variation. It is not uncertainty about whether the assumed trial model is scientifically adequate, nor uncertainty in a model parameter. It answers a narrower question: if the same simulation procedure were repeated with a new set of random draws, how much might its estimated success probability vary? Quantifying this error helps determine whether a reported probability is precise enough for its intended decision.

Definition: Let \(p\) be the probability that one replicate meets a specified success criterion under a fixed scenario and simulation procedure. For \(R\) independent replicates, define \(I_r=1\) if replicate \(r\) succeeds and \(I_r=0\) otherwise. The number of successes is \(K=\sum_{r=1}^{R} I_r\), and the estimated success probability is \(\hat{p}=K/R\). The Monte Carlo standard error (MCSE) estimates the standard deviation of \(\hat{p}\) across repeated sets of \(R\) replicates.

Monte Carlo standard error

For a fixed scenario and success rule, each \(I_r\) has a Bernoulli distribution with probability \(p\). If replicate outcomes are independent and generated under the same procedure, \(K\) has a binomial distribution with parameters \(R\) and \(p\). Thus \(E(\hat{p})=p\), and

\operatorname{Var}(\hat{p})=\operatorname{Var}\left(\frac{K}{R}\right) =\frac{p(1-p)}{R}.

Since \(p\) is unknown, estimate it with \(\hat{p}\). The usual estimated MCSE is

\widehat{\operatorname{SE}}_{\mathrm{MC}}(\hat{p}) =\sqrt{\frac{\hat{p}(1-\hat{p})}{R}}.

The estimate is on the probability scale. Multiplying by 100 expresses it in percentage points. It is not the standard deviation of individual participant outcomes; it is the estimated standard deviation of the replicate-based probability estimate.

Worked Example: Calculating an MCSE

Suppose a hypothetical trial-simulation scenario produces 72 successful replicates out of \(R=1{,}000\). The estimated success probability is

\hat{p}=\frac{K}{R}=\frac{72}{1{,}000}=0.072.

The estimated MCSE is

\widehat{\operatorname{SE}}_{\mathrm{MC}} =\sqrt{\frac{0.072(1-0.072)}{1{,}000}} =\sqrt{0.000066816} =0.008174.

So the estimated MCSE is 0.008174 on the probability scale, or 0.8174 percentage points. This quantifies the typical sampling variation from finite simulation replicates; it does not mean that 0.8174 percentage points is a bound on the error in this particular estimate.

Confidence intervals for a simulated probability

An approximate 95% interval can be formed as \(\hat{p}\pm1.96\widehat{\operatorname{SE}}_{\mathrm{MC}}\). This Wald interval is simple, but can perform poorly when the replicate count is small or the estimated probability is near 0 or 1. It can even extend beyond the possible range from 0 to 1. A Wilson interval is generally preferable for a binomial probability.

To obtain the Wilson interval, start with the standardized binomial-proportion statistic and retain values of \(p\) that satisfy the two-sided 95% acceptance condition. Inverting that condition gives the Wilson limits below, where \(z=1.96\) for a 95% interval. The center and half-width form makes clear that both depend on the observed proportion and replicate count.

\text{center} =\frac{\hat{p}+\frac{z^2}{2R}}{1+\frac{z^2}{R}}, \qquad \text{half-width} =\frac{z}{1+\frac{z^2}{R}} \sqrt{\frac{\hat{p}(1-\hat{p})}{R}+\frac{z^2}{4R^2}}.

\text{Wilson interval} =\text{center}\ \pm\ \text{half-width}.

The interval procedure has approximately 95% coverage under repeated independent binomial simulation for a fixed \(p\). In context, it describes uncertainty from the finite number of simulation replicates. It does not assign a 95% probability that the fixed probability \(p\) lies in the particular interval observed.

Worked Example: A Wilson Interval for Success Probability

In a hypothetical simulation, \(K=82\) of \(R=100\) replicates meet the success criterion, so \(\hat{p}=0.82\). For a 95% Wilson interval, use \(z=1.96\), with \(z^2=3.8416\). The denominator is

1+\frac{z^2}{R}=1+\frac{3.8416}{100}=1.038416.

The center is

\frac{0.82+\frac{3.8416}{200}}{1.038416} =\frac{0.839208}{1.038416} =0.8082.

The half-width is

\frac{1.96}{1.038416} \sqrt{\frac{0.82(0.18)}{100}+\frac{3.8416}{4(100)^2}} = \frac{1.96}{1.038416}\sqrt{0.00157204} =0.07484.

Therefore, the 95% Wilson interval is approximately \(0.8082\pm0.07484\), or \((0.7333,\,0.8830)\). The simulated success probability is estimated as 82%; under repeated simulation with the same fixed setup, the Wilson procedure has approximately 95% coverage for the underlying success probability.

Choosing the number of replicates

Replicate counts can be planned around a target precision, often expressed as a desired confidence-interval half-width \(d\). The approximate normal interval has half-width \(z\sqrt{p(1-p)/R}\). Solving for \(R\) gives a planning calculation. Because \(p\) is not known before simulation, a preliminary estimate can be used, or the conservative bound \(p(1-p)\leq0.25\) can be used when a guarantee against the worst-case binomial variance is preferred.

d\approx z\sqrt{\frac{p(1-p)}{R}} \quad\Longrightarrow\quad R\approx\frac{z^2p(1-p)}{d^2}.

p(1-p)\leq0.25 \quad\Longrightarrow\quad R\geq\frac{z^2}{4d^2} \quad\text{for the conservative normal-approximation plan.}

These are planning approximations, not promises that the realized interval will have exactly the requested half-width. The actual interval width depends on the observed proportion. A practical plan specifies a target, uses a preliminary or conservative probability for the calculation, rounds \(R\) up to a whole number, and then checks the interval width after simulation.

Worked Example: Planning Replicates for a Target Half-Width

Suppose a preliminary run suggests a hypothetical success probability near 0.80. For a 95% interval with a target approximate half-width of 0.01 (one percentage point), the plug-in planning calculation is

R\approx\frac{1.96^2(0.80)(0.20)}{0.01^2} =\frac{3.8416(0.16)}{0.0001} =6{,}146.56.

Round up to \(R=6{,}147\) replicates. This is an approximate plan based on the preliminary probability, not a guarantee that the final interval will have a half-width no greater than 0.01. If the probability were unknown and the conservative variance bound were used instead, the calculation would be

R\geq\frac{1.96^2}{4(0.01)^2} =\frac{3.8416}{0.0004} =9{,}604.

The larger count protects against the worst-case binomial variance under this approximation. After running the chosen number of replicates, compute the interval from the observed \(K/R\) and verify whether its half-width meets the precision requirement.

Assumptions, diagnostics, and alternatives

The binomial calculations apply when each replicate has a clearly defined binary success outcome, all replicates target the same fixed scenario and decision rule, and replicate outcomes are independent with the same success probability. If the procedure changes between replicates, or if outcomes are dependent, the binomial variance and interval formulas may not be appropriate. In particular, the independence assumption concerns the replicate outcomes used to estimate one probability.

Useful checks include confirming that \(0\leq K\leq R\), that the reported estimate equals \(K/R\), and that the interval is within \([0,1]\). Inspect the estimated probability and interval width together: a small MCSE in absolute terms may still be inadequate for a rare success probability, while a probability near a boundary makes a Wald interval especially questionable. If a probability is very rare, also consider whether the run produced enough successes to support a stable estimate; an observed zero does not establish that the underlying probability is zero.

The MCSE formula is an estimate; substituting \(\hat{p}=0\) or \(\hat{p}=1\) yields zero, which can misleadingly suggest certainty after observing all failures or all successes. A Wilson interval avoids that degeneracy and remains within the probability range. An exact binomial interval is another option, particularly when counts are very small; it is obtained by inverting binomial tail tests and is generally conservative. The Wald interval is easiest to calculate, but its weak behavior near the boundaries is a reason not to rely on it as the default.

Increasing \(R\) reduces MCSE approximately in proportion to \(1/\sqrt{R}\), not \(1/R\). Thus, cutting MCSE in half requires about four times as many replicates. More replicates reduce Monte Carlo error but do not repair a poor scenario definition, an inappropriate success criterion, or an unrealistic simulation model. Those are distinct sources of uncertainty and should not be represented by the binomial interval.

Common pitfall: Treating the observed success fraction as exact, or interpreting its MCSE as uncertainty in the scientific model. The correction is to report the replicate count, estimate, and an interval or MCSE for finite-simulation variation, while keeping that variation distinct from uncertainty about the assumptions used to generate the trials.

Check Your Understanding

Work each problem, then compare your reasoning with the solutions.

A hypothetical simulation has 150 successes in 2,500 independent replicates. Calculate \(\hat{p}\) and its estimated MCSE.
Why is a Wald interval a poor choice when the estimated success probability is 0 or very close to 0?
Using a 95% normal approximation and preliminary \(p=0.50\), calculate the replicate count for an approximate half-width of 0.02.
Does an MCSE quantify uncertainty about whether the trial model represents clinical reality? Explain.
Why can observing zero successes in a finite simulation not establish that the underlying success probability is zero?

Solutions. (1) \(\hat{p}=150/2{,}500=0.0600\). The MCSE is \(\sqrt{0.06(0.94)/2{,}500}=\sqrt{0.00002256}=0.004750\), or 0.4750 percentage points. (2) The Wald interval is centered at the observed proportion with a symmetric normal-based margin; near zero it can extend below zero, and at zero its plug-in standard error is zero. A Wilson or exact binomial interval is more appropriate. (3) \(R\approx1.96^2(0.50)(0.50)/0.02^2=3.8416(0.25)/0.0004=2{,}401\) replicates. (4) No. MCSE describes variation in the estimated probability due to a finite number of simulation replicates under a fixed simulation setup. It does not measure model or clinical-relevance uncertainty. (5) A finite set of Bernoulli draws can contain no successes even when \(p>0\); zero observed successes gives \(\hat{p}=0\), not proof that \(p=0\).