Separating Sampling and Simulation Uncertainty

Why Separate the Two Sources of Variation?

A simulated estimator depends on the observed sample and on numerical randomness used to approximate an expectation, probability, or objective function. Repeating the study with a new sample changes the data; rerunning the algorithm with the same data changes the simulation draws. These are distinct sources of uncertainty. If they are combined without being identified, it can be difficult to tell whether an imprecise estimate reflects limited data or avoidable Monte Carlo noise.

This distinction matters for inference. A standard error intended to describe repeated-sample variability should not omit substantial simulation variability. Conversely, simulation variation should not be mistaken for variation that would remain even if the simulation were made arbitrarily precise. The decomposition below gives a precise way to describe both contributions and to study them with repeated datasets and repeated runs.

Definition: Let \(D\) denote the random dataset and \(U\) the random simulation draws used by an estimator \(\hat{\theta}(D,U)\) for a parameter \(\theta\). Its simulation component is \(E_D[\operatorname{Var}_U(\hat{\theta}\mid D)]\). Its data component is \(\operatorname{Var}_D(E_U[\hat{\theta}\mid D])\). The data component describes variation across datasets in the estimator's conditional simulation mean; it equals the variance of the corresponding exact, non-simulated estimator when that conditional mean equals that estimator.

The distinction in the last sentence is important. A simulation procedure can have conditional bias: for a particular dataset, its expected value over simulation draws may differ from the estimator that would be obtained without simulation. In that case, the data component defined above is the variance of the simulated estimator's conditional mean, not necessarily the sampling variance of an exact estimator. Simulation bias is a separate issue from simulation variance.

The Variance Decomposition

The law of total variance gives an exact identity for any estimator with finite variance:

\operatorname{Var}(\hat{\theta}) = \operatorname{Var}_D\!\left(E_U[\hat{\theta}\mid D]\right) + E_D\!\left[\operatorname{Var}_U(\hat{\theta}\mid D)\right].

The first term is variation in the estimator's conditional mean as the dataset changes. The second is the average variation across simulation runs while holding the dataset fixed. This identity does not require the estimator to be unbiased. It describes variance around \(E[\hat{\theta}]\); it does not by itself describe mean squared error about \(\theta\).

Deriving the Decomposition

Write \(m(D)=E_U[\hat{\theta}\mid D]\). Add and subtract this conditional mean from the centered estimator:

\hat{\theta}-E[\hat{\theta}] = \{\hat{\theta}-m(D)\}+\{m(D)-E[\hat{\theta}]\}.

The first bracket has conditional expectation zero given \(D\). Consequently its covariance with any function of \(D\), including the second bracket, is zero:

E\!\left[(\hat{\theta}-m(D))(m(D)-E[\hat{\theta}])\right] = E\!\left[E[\hat{\theta}-m(D)\mid D](m(D)-E[\hat{\theta}])\right] =0.

Squaring the centered expression, taking expectations, and using this zero cross-product gives

\operatorname{Var}(\hat{\theta}) = E[(\hat{\theta}-m(D))^2] + E[(m(D)-E[\hat{\theta}])^2] = E_D[\operatorname{Var}_U(\hat{\theta}\mid D)] + \operatorname{Var}_D(m(D)).

This derivation also explains why the simulation variance is averaged over datasets: simulation variability may differ from one dataset to another.

Worked Example: Combining Components for an Interval

Suppose a hypothetical estimator of a treatment-effect parameter, measured in percentage points, has data variance \(0.64\) and simulation variance \(0.20\). Assume the conditional simulation mean agrees with the exact estimator, the estimator has negligible bias relative to its standard error, and the standardized error about \(\theta\) is approximately standard normal. Find the total standard error and a 95% interval when \(\hat{\theta}=2.40\).

The variance decomposition gives total variance \(0.64+0.20=0.84\), so

SE=\sqrt{0.84}=0.9165\ \text{percentage points}.

Using the standard normal critical value \(1.96\), the interval is

2.40\pm 1.96(0.9165) =2.40\pm1.796 =[0.604,\ 4.196]\ \text{percentage points}.

Under the stated centering and normal-approximation conditions, this procedure has approximately 95% repeated-sampling coverage when both data and simulation randomness are represented. If simulation variance were omitted, the reported standard error would be \(\sqrt{0.64}=0.800\). Its interval would be too narrow for the stated total variance: under a centered normal error with variance \(0.84\), its coverage would be about \(2\Phi(1.96(0.800)/0.9165)-1=91.3\%\).

Normality and a correct variance estimate alone are not enough for the 95% coverage claim. For example, if the estimator's error were distributed as \(N(1,0.84)\), its standard deviation would still be \(\sqrt{0.84}\), but the interval centered at \(\hat{\theta}\) would miss the target much more often because the error is biased. The centering condition must be justified as well.

Estimating Components with Repeated Runs

A nested repeated-run design uses \(M\) independently generated datasets and \(R\) independent simulation runs for each dataset. Denote the result from dataset \(m\) and run \(r\) by \(\hat{\theta}_{mr}\). For each dataset, calculate its average over simulation runs, \(\bar{\theta}_m=R^{-1}\sum_{r=1}^R\hat{\theta}_{mr}\), and the grand average \(\bar{\theta}_{..}=M^{-1}\sum_{m=1}^M\bar{\theta}_m\).

When the design is balanced, simulation runs are independent conditional on each dataset, and the conditional simulation variance is common across datasets, the within-dataset and between-dataset mean squares estimate the two components. Define

MS_W= \frac{\sum_{m=1}^M\sum_{r=1}^R(\hat{\theta}_{mr}-\bar{\theta}_m)^2}{M(R-1)}, \qquad MS_B= \frac{R\sum_{m=1}^M(\bar{\theta}_m-\bar{\theta}_{..})^2}{M-1}.

The expected within mean square is the simulation variance, \(\sigma^2_{\mathrm{sim}}\). The expected between mean square is \(\sigma^2_{\mathrm{sim}}+R\sigma^2_{\mathrm{data}}\), because each dataset average still contains simulation variance reduced by a factor of \(R\). Thus the method-of-moments estimates are

\widehat{\sigma^2_{\mathrm{sim}}}=MS_W, \qquad \widehat{\sigma^2_{\mathrm{data}}}=\frac{MS_B-MS_W}{R}.

In a small or unbalanced study these estimates can be noisy; the data-component estimate can even be negative. A negative estimate is a warning that the design provides little evidence for a positive between-dataset component, not a claim that a true variance is negative. The common-variance and independence conditions also matter: if simulation variance differs considerably by dataset, or runs share random draws in a way that induces dependence, this simple calculation need not estimate the intended components.

Worked Example: A Small Nested Simulation Study

Consider a hypothetical experiment with \(M=4\) independently generated datasets and \(R=2\) independent simulation runs per dataset. The estimated effects, in score points, are:

Dataset	Run 1	Run 2	Dataset average
1	1.0	1.2	1.1
2	1.8	2.0	1.9
3	2.2	2.4	2.3
4	2.8	3.0	2.9

The grand average is \(\bar{\theta}_{..}=(1.1+1.9+2.3+2.9)/4=2.05\). Each dataset has two deviations from its average, \(-0.1\) and \(0.1\), so the total within sum of squares is \(4(0.1^2+0.1^2)=0.08\). Therefore,

MS_W=\frac{0.08}{4(2-1)}=0.0200\ \text{score points}^2.

The deviations of the dataset averages from the grand average are \(-0.95,-0.15,0.25,0.85\). Their squared sum is \(0.9025+0.0225+0.0625+0.7225=1.71\). Hence

MS_B=\frac{2(1.71)}{4-1}=1.14\ \text{score points}^2, \qquad \widehat{\sigma^2_{\mathrm{data}}}=\frac{1.14-0.0200}{2}=0.560.

The estimated variance components are \(0.560\) score points squared for the data component and \(0.0200\) score points squared for simulation. For one simulation run per dataset, the estimated total variance is \(0.580\) score points squared, with standard deviation \(\sqrt{0.580}=0.7616\) score points. The estimate of the simulation share is \(0.0200/0.580=3.45\%\). This small illustrative design suggests relatively little simulation variability, but four datasets provide limited precision; the result should not be treated as a definitive variance estimate.

Deciding Whether Simulation Error Matters

A variance share is useful, but it is not the only relevant measure. Simulation noise may be small compared with total variability yet still matter for a decision near a significance threshold, or when comparing two close estimates. Assess its scale in the units of the parameter, its contribution to the standard error, and the intended inferential use. Also examine whether increasing the number of draws reduces the simulation component as expected. For many Monte Carlo averages based on independent draws, simulation variance is approximately proportional to \(1/K\), where \(K\) is the number of draws; changing the design or using dependent draws can alter that rate or its constant.

Worked Example: Assessing the Value of More Draws

Suppose a hypothetical estimator has data variance \(0.36\) and simulation variance \(0.09\) when it uses \(K=100\) independent draws. Assume the simulation variance scales as \(1/K\), with the data component unchanged. Compare the total standard errors at \(K=100\) and \(K=900\).

At \(K=100\), the total variance is \(0.36+0.09=0.45\), and the standard error is \(\sqrt{0.45}=0.6708\) units of the parameter. Increasing the draw count ninefold predicts simulation variance \(0.09(100/900)=0.0100\). At \(K=900\), total variance is \(0.36+0.0100=0.370\), and the standard error is \(\sqrt{0.370}=0.6083\) units.

The simulation variance share falls from \(0.09/0.45=20.0\%\) to \(0.0100/0.370=2.70\%\). The standard error falls by about \(0.0625\) units, or \(9.32\%\) relative to its value at \(K=100\). More simulation draws reduce simulation noise, but they cannot remove the data component. Whether the reduction is consequential depends on the precision needed for the analysis, not on a universal cutoff.

From Variance to Inference

For inference about \(\theta\), the total variance is not the whole story. The mean squared error satisfies

E[(\hat{\theta}-\theta)^2] = \operatorname{Var}(\hat{\theta})+\{E[\hat{\theta}]-\theta\}^2.

Thus a variance-based interval centered at \(\hat{\theta}\) needs the estimator to be approximately centered at \(\theta\), or to have bias negligible relative to its standard error. A commonly used sufficient description is that the standardized error, including both sources of randomness, is approximately standard normal:

\frac{\hat{\theta}-\theta}{\sqrt{\sigma^2_{\mathrm{data}}+\sigma^2_{\mathrm{sim}}}} \ \approx\ N(0,1).

When this approximation is appropriate and the component estimates are reliable, a corresponding approximate 95% interval is \(\hat{\theta}\pm1.96\sqrt{\widehat{\sigma^2_{\mathrm{data}}}+\widehat{\sigma^2_{\mathrm{sim}}}}\). If the data and simulation components are estimated from a finite nested experiment, their estimation uncertainty may also matter; the displayed interval treats the estimated variance as sufficiently accurate.

A frequent pitfall is to claim that a normal estimator and a correct standard error automatically give 95% coverage. They do not: if the estimator is biased relative to the target by an amount that is not negligible compared with its standard error, a nominal interval can undercover. Another mistake is to increase the simulation count and assume that all uncertainty must disappear. More draws address simulation variability, not sampling variation, and they do not automatically correct simulation bias.

Key takeaway: Separate the estimator's variation across datasets from its variation across simulation runs on the same dataset. Use the total variance for inference only with justified centering, an appropriate distributional approximation, and variance estimates that represent both sources of randomness.

Check Your Understanding

Use the definitions and conditions in this tutorial to answer the following.

In the variance decomposition, what is held fixed when calculating the simulation component?
Why does the between-dataset mean square contain simulation variance when each dataset average is based on multiple runs?
A simulation-aware standard error is estimated correctly, but the estimator has bias that is large relative to that standard error. What implication does this have for a nominal 95% interval?
If simulation variance is one quarter of total variance for a single-run estimate, what proportion of total variance is the data component?
When is the data component from the conditional-variance decomposition equal to the sampling variance of an exact, non-simulated estimator?