Why Compare Simulation Estimates Carefully?
A simulation study often compares two hypothetical alternatives: for example, two deductible designs, two claim-handling rules, or two proposed operating policies. The quantity of interest is usually the difference between their expected outputs. Because each output is estimated with simulation, that difference is also subject to Monte Carlo error. A difference that appears positive in one run may be small relative to the uncertainty in the estimates.
This tutorial considers independent simulations: the replications used for one alternative are independent of those used for the other. The objective is to quantify simulation uncertainty in the difference, not uncertainty about whether the assumed model itself is appropriate. The resulting interval does not include parameter-estimation error, model risk, or uncertainty in inputs held fixed during the simulation.
The order in the subtraction matters. If \(D\) is defined as alternative A minus alternative B, then a positive estimate means that A has the higher average output. Whether that is favorable depends on what the output measures: a higher profit may be desirable, while a higher claim cost may not be.
Standard Error for an Independent Difference
For \(n_A\) replications of A and \(n_B\) replications of B, let \(s_A^2\) and \(s_B^2\) be the usual sample variances of their outputs. The estimated standard error of \(D\) combines the uncertainty from both sample means.
The plus sign is essential. The two estimates may be subtracted to form the difference, but independent estimation errors contribute variance from both sides. The formula allows unequal sample sizes and unequal output variances; it does not require a common variance.
Calculate \(D=\bar{X}_A-\bar{X}_B\), keeping the chosen order consistent with the interpretation.
Calculate \(s_A^2/n_A\) and \(s_B^2/n_B\), add them, and take the square root.
For sufficiently large samples, a common approximate \(100(1-\alpha)\%\) interval is \(D\pm z_{1-\alpha/2}\operatorname{SE}(D)\), where \(z_{1-\alpha/2}\) is the corresponding standard normal quantile.
State the interval in the output's units and explain what values of the true difference it supports.
Why the Variances Add
A brief derivation shows why the standard error has this form. By linearity of expectation, \(E[D]=E[\bar{X}_A]-E[\bar{X}_B]=\mu_A-\mu_B=\Delta\), so the difference of sample means estimates the desired difference without bias. For any two random variables, the variance of a difference is the sum of their variances minus twice their covariance. Independence makes that covariance zero.
The last equality uses independence between alternatives and the fact that the variance of a sample mean of independent, identically distributed outputs is the output variance divided by the replication count. Replacing each unknown variance with its sample estimate gives the estimated standard error. If the two sets of replications are not independent, the covariance term need not be zero, so this independent-simulation formula does not apply.
Confidence Intervals and Their Assumptions
When both replication counts are reasonably large and the outputs have finite variances, the central limit theorem often makes the difference of sample means approximately normal. A 95% interval then uses \(z_{0.975}\approx1.96\). The interval describes the Monte Carlo precision of the estimate: under repeated independent simulation runs, intervals constructed this way will capture the fixed target difference about 95% of the time, subject to the approximation.
For smaller samples, or when a normal approximation is not credible, a Welch interval is often preferable if the replication outputs are approximately normal. It uses the same estimated standard error but a \(t\) critical value and an estimated degrees of freedom:
The Welch interval is \(D\pm t_{1-\alpha/2,\nu}\operatorname{SE}(D)\). Its degrees of freedom need not be an integer; software or a \(t\)-table can be used to obtain the critical value. It avoids assuming the two alternatives have equal output variances. Neither the normal nor Welch interval automatically works well for very small samples from highly skewed or heavy-tailed output distributions.
Worked Example: Comparing Average Claim Costs
A hypothetical insurer independently simulates average annual claim cost under two coverage designs. Design A has \(n_A=400\), sample mean \(\bar{X}_A=\$1{,}218\), and sample standard deviation \(s_A=\$120\). Design B has \(n_B=625\), sample mean \(\bar{X}_B=\$1{,}200\), and sample standard deviation \(s_B=\$100\). Estimate the difference A minus B and form an approximate 95% confidence interval. Treat the replication counts as large enough for a normal approximation.
Method. Use the independent-sample standard error and a normal confidence interval. The estimated difference is \(\$1{,}218-\$1{,}200=\$18\).
The two estimated variance contributions are \(s_A^2/n_A=(\$120)^2/400=\$14{,}400/400=\$36\) and \(s_B^2/n_B=(\$100)^2/625=\$10{,}000/625=\$16\). Therefore,
Using \(z_{0.975}=1.96\), the margin of error is \(1.96(\$7.2111)=\$14.1338\). Thus the interval is \(\$18\pm\$14.1338\), or \((\$3.87,\$32.13)\) after rounding to the nearest cent. This interval is entirely positive, indicating that the simulation estimates support a higher mean claim cost for A than B, while the interval width shows the remaining Monte Carlo uncertainty.
Check. The midpoint is \((\$3.87+\$32.13)/2=\$18.00\), which agrees with the estimated difference. The method also avoids a common mistake: treating the difference as if only one alternative had simulation error.
Binary Outputs: Comparing Simulated Probabilities
A simulated event indicator is a numerical output equal to 1 when the event occurs and 0 otherwise. Its sample mean is the estimated event probability. For independent Bernoulli replications, the estimated variance of a sample mean can be written \(\hat{p}(1-\hat{p})/n\), where \(\hat{p}\) is the observed event proportion. The independent-difference formula therefore applies directly to two estimated probabilities. A normal interval may be inaccurate when event counts are small or estimated probabilities are near 0 or 1; in such cases, a suitable interval method for proportions should be considered instead.
Worked Example: Comparing Two Simulated Event Probabilities
In a hypothetical risk analysis, 80 of 1,000 independent simulations under Plan A produce a specified loss event, while 90 of 1,500 independent simulations under Plan B produce it. Estimate \(P(A\text{ event})-P(B\text{ event})\) and construct an approximate 95% interval using the normal approximation.
Method. Estimate the two probabilities by their observed proportions: \(\hat{p}_A=80/1{,}000=0.0800\) and \(\hat{p}_B=90/1{,}500=0.0600\). The estimated difference is \(0.0800-0.0600=0.0200\).
The estimated variance contributions are \(0.0800(1-0.0800)/1{,}000=0.0736/1{,}000=0.0000736\) and \(0.0600(1-0.0600)/1{,}500=0.0564/1{,}500=0.0000376\). Hence the standard error is \(\sqrt{0.0000736+0.0000376}=\sqrt{0.0001112}=0.010545\).
The 95% margin of error is \(1.96(0.010545)=0.020668\). The interval is \(0.0200\pm0.020668=(-0.000668,0.040668)\), or \((-0.0007,0.0407)\) to four decimal places. The interval includes zero, so this simulation does not establish a positive difference at this confidence level.
Check. The interval is centered at \(0.0200\), and its endpoints differ from the center by approximately \(0.0207\). The method avoids confusing the estimated event counts with the estimated probabilities: each count must be divided by its own replication count.
Welch Interval with Unequal Variances
The large-sample normal interval is convenient, but a small-sample comparison can call for a \(t\) critical value. The following example uses unequal sample sizes and sample variances, so it does not pool the two variance estimates.
Worked Example: Comparing Mean Processing Times
A hypothetical analyst independently simulates processing time, measured in thousands of dollars of system cost, for two claim workflows. For A, \(n_A=16\), \(\bar{X}_A=52\), and \(s_A=4\). For B, \(n_B=25\), \(\bar{X}_B=49\), and \(s_B=3\). Assume the replication outputs are approximately normal. Construct a 95% Welch interval for A minus B.
Method. First, \(D=52-49=3\). The variance contributions are \(s_A^2/n_A=16/16=1\) and \(s_B^2/n_B=9/25=0.36\). Thus \(\operatorname{SE}(D)=\sqrt{1+0.36}=\sqrt{1.36}=1.1662\).
The Welch degrees of freedom are \(\nu=(1.36)^2/\{1^2/15+0.36^2/24\}=1.8496/(0.066667+0.005400)=25.66\), approximately 26. A 95% \(t\) critical value for about 26 degrees of freedom is approximately 2.06. The margin is \(2.06(1.1662)=2.4024\), so the interval is \(3\pm2.4024=(0.598,5.402)\) in thousands of dollars, rounded to three decimal places.
Check. The interval is centered at 3 and lies above zero. It uses the larger uncertainty contribution from each alternative separately; it does not require the alternatives to share a common variance. The method avoids pooling the variances merely because the two means are being compared.
Common Mistakes and Interpretation
A frequent error is to subtract the variance contributions, or to treat the standard error of the difference as the standard error from one alternative alone. Independent estimation errors add in variance, even though the estimated outputs are subtracted. For the claim-cost example, the incorrect calculation \(\sqrt{(14{,}400+10{,}000)/(400+625)}=\$4.88\) does not give the standard error of the difference: it combines the variance estimates using a pooled-style denominator that is not appropriate for this difference. The correct calculation is \(\sqrt{14{,}400/400+10{,}000/625}=\$7.21\).
Also check the interpretation of zero. If a confidence interval for \(\Delta\) contains zero, the data do not clearly distinguish the alternatives at the stated confidence level; this does not prove their true means are identical. If the interval excludes zero, it gives evidence about the direction of the mean difference under the simulation assumptions, but it does not by itself establish that the size of the difference is practically important.
Finally, the formula in this tutorial relies on independent simulation runs across alternatives. If randomness is shared between alternatives, the covariance is generally not zero, and the independent-sample standard error is not the appropriate calculation. The output should also be a sample mean, or a quantity for which the stated approximation has been justified. A confidence interval for a simulated quantile, for example, requires a method suited to quantile uncertainty rather than simply applying this mean formula.
Check Your Understanding
For each problem, define the difference as alternative A minus alternative B. Unless a problem says otherwise, use the normal approximation and report interval endpoints to three decimals.
- Two independent simulations estimate mean repair costs. A has \(n_A=400\), mean \(\$14.2\), and standard deviation \(\$8\); B has \(n_B=625\), mean \(\$13.5\), and standard deviation \(\$10\). Find the standard error and a 90% interval using \(z=1.645\). (Answer: \(D=\$0.700\), \(\operatorname{SE}=\$0.566\), 90% interval \((-\$0.231,\$1.631)\).)
- In independent simulations of a hypothetical event, A has 60 events in 500 replications and B has 72 events in 800 replications. Find the estimated probability difference and a 95% Wald interval. (Answer: \(D=0.0300\), \(\operatorname{SE}\approx0.0177\), 95% interval \((-0.0047,0.0647)\).)
- For independent simulated outputs, A has mean 76, standard deviation 15, and \(n_A=100\); B has mean 71, standard deviation 12, and \(n_B=144\). Find the standard error and approximate 95% interval using \(z=1.96\). (Answer: \(D=5.000\), \(\operatorname{SE}=1.803\), 95% interval \((1.467,8.533)\).)
- Two independent simulations estimate mean daily costs. A has mean 3.4, standard deviation 0.8, and \(n_A=40\); B has mean 3.1, standard deviation 0.6, and \(n_B=60\). Find a 90% interval using \(z=1.645\), and state whether it includes zero. (Answer: \(D=0.300\), \(\operatorname{SE}=0.148\), 90% interval \((0.056,0.544)\); it does not include zero.)
- Independent simulations give A a mean output of 22, sample variance 9, and \(n_A=36\), and B a mean output of 19, sample variance 16, and \(n_B=64\). Find the estimated standard error and approximate 95% interval using \(z=1.96\). (Answer: \(D=3.000\), \(\operatorname{SE}=0.707\), 95% interval \((1.614,4.386)\).)