Why Verify Against Known Results?
A simulation can produce plausible-looking operating characteristics while implementing the wrong data-generating model, decision rule, or summary calculation. Verification asks a focused question: when a component has a result that can be calculated independently, does the code reproduce that result? The benchmark may be an exact probability, a hand-calculated expectation, a known variance, or a mathematically justified limit.
This is different from estimating the performance of a full adaptive design, where the answer may not be available in closed form. Verification uses simpler cases to establish that important pieces behave as specified before they are combined. Unit tests of adaptive decision rules check particular inputs and branches; analytical benchmarks add an external numerical target. Neither alone proves that a complete design is correct.
A comparison must account for two sources of discrepancy. Monte Carlo variation is the random difference between a finite simulation and its exact target. An implementation error is a systematic or structural mismatch, such as using the wrong success probability or applying a stopping condition at the wrong time. A single simulated estimate need not equal its target exactly.
Choose Benchmarks That Match the Component
A useful benchmark is simple enough to calculate reliably and specific enough to expose a plausible coding error. For a generator, check its distribution or moments. For an estimator, check its expected value and variance under a tractable model. For an adaptive rule, enumerate possible interim outcomes and calculate the probability of each resulting decision or sample size. For a limiting check, state what quantity approaches what value, and under which assumptions.
Start by stating the model. For example, if observations \(X_1,\ldots,X_n\) are independent with means \(\mu_i\) and finite variances \(\sigma_i^2\), then the mean \(\bar X=(X_1+\cdots+X_n)/n\) has expectation equal to the average of the observation means. Independence is needed to add variances without covariance terms; equal variances are not required.
The variance result follows by applying the variance rules to the weighted sum. Since the observations are independent, every covariance between distinct observations is zero; scaling each \(X_i\) by \(1/n\) scales its variance by \(1/n^2\). If, in addition, all observations have the same variance \(\sigma^2\), the expression reduces to \(\operatorname{Var}(\bar X)=\sigma^2/n\). That simplified formula must not be used merely because observations are independent.
This also gives a careful limiting check. If all variances equal \(\sigma^2<\infty\), then \(\sigma^2/n\) tends to zero as \(n\) grows. More generally, for independent observations with possibly unequal variances, the variance tends to zero when \(\sum_{i=1}^n\sigma_i^2/n^2\) tends to zero. For example, a uniform upper bound on the variances is sufficient: if \(\sigma_i^2\le C\) for every \(i\), then the variance is at most \(C/n\). Independence and finite variances alone, without a condition on how those variances grow, do not guarantee this limit.
Worked Example: Checking a Mean With Unequal Variances
Suppose a hypothetical simulation generates two independent observations, \(X_1\sim N(10,4)\) and \(X_2\sim N(12,16)\), and reports their average. The second argument of each normal distribution is its variance. Find the exact mean and variance of the reported average, then compare them with a hypothetical simulation output.
The model specifies independent observations, with finite variances 4 and 16. Their average has expected value \(E(\bar X)=(10+12)/2=11\). Its variance is \((4+16)/2^2=20/4=5\), so its standard deviation is \(\sqrt{5}=2.236\).
Now suppose the code repeats this experiment \(R=100{,}000\) times and reports an average of 10.994 across the simulated averages. The Monte Carlo standard error of that reported average is \(\sqrt{5/100{,}000}=0.007071\), because the simulated average in each replication has variance 5. The observed difference from the benchmark is \(10.994-11=-0.006\), or about \(-0.85\) Monte Carlo standard errors. This is consistent with ordinary simulation variation; it is not evidence by itself of a bug.
Using the equal-variance shortcut would give the wrong target here. There is no single common \(\sigma^2\) for these two observations; the correct variance is 5, not a formula based on treating them as identically distributed.
Comparing a Simulated Rate With an Exact Probability
For a binary outcome, suppose each replication produces an indicator \(I_r\) that equals 1 when an event occurs and 0 otherwise. Under the specified model, if the event probability is \(q\), the \(R\) independent replications give a simulated rate \(\hat q=R^{-1}\sum_{r=1}^R I_r\). Its expectation is \(q\), and its variance is \(q(1-q)/R\). Thus its Monte Carlo standard error is \(\sqrt{q(1-q)/R}\). This formula assumes independent replications with the same event probability; dependence between replications requires a different calculation.
When \(q\) is not known, the usual estimated Monte Carlo standard error substitutes \(\hat q\) for \(q\). When an exact benchmark is available, using \(q\) directly is useful for evaluating the expected simulation variation around that target. A discrepancy of a few Monte Carlo standard errors is a diagnostic signal, not an automatic proof of an error: the calculation is approximate for finite \(R\), and repeated checks can produce occasional large discrepancies by chance. A gross or repeatable mismatch calls for investigation.
Worked Example: Verifying a Bernoulli Event Probability
Consider hypothetical code that generates four independent Bernoulli outcomes, each with success probability \(p=0.25\), and records whether there are at least two successes. Calculate the exact event probability and use it to assess a simulated rate from \(R=40{,}000\) replications.
The number of successes has a binomial distribution with four trials and success probability 0.25. The event is the complement of zero or one success: \(q=1-P(0)-P(1)=1-(0.75)^4-4(0.25)(0.75)^3=0.26171875\).
Suppose the code records the event in 10,500 of its 40,000 replications. Its simulated rate is \(10{,}500/40{,}000=0.2625\). Under independent replications, the benchmark Monte Carlo standard error is \(\sqrt{0.26171875(1-0.26171875)/40{,}000}=0.002198\). The difference is \(0.2625-0.26171875=0.0007813\), about 0.36 standard errors. That difference is small relative to expected Monte Carlo variation, so this result is consistent with the exact benchmark.
This check can detect, for example, a generator using 0.2 instead of 0.25 or a condition that mistakenly counts at least three successes. It does not establish that every aspect of the generator is correct; additional benchmarks should target other relevant features.
Verify an Adaptive Rule by Enumerating Its Paths
For a small adaptive design, exact enumeration can provide a direct benchmark for both a decision probability and a resource-use summary. The calculation requires a specified outcome model and a precise rule: when the interim occurs, what triggers stopping or adaptation, and how final success is defined. For binary outcomes, a small number of patients can make every possible path enumerable by hand.
Worked Example: An Interim Futility Rule
Consider a hypothetical two-stage trial with independent Bernoulli outcomes, each having success probability \(p=0.40\). The first stage enrolls two patients. If there are zero successes, the trial stops for futility. Otherwise, two more patients are enrolled. The trial is ultimately called successful if at least three of the four possible patients would have succeeded. Find the exact probability of trial success and the expected sample size.
Write \(K_1\) for the number of successes among the first two patients. Success is impossible if \(K_1=0\), because the trial stops. If \(K_1=1\), both later patients must succeed, which has probability \((0.40)^2=0.16\). The probability of \(K_1=1\) is \(2(0.40)(0.60)=0.48\), so this path contributes \(0.48(0.16)=0.0768\).
If \(K_1=2\), at least one of the two later patients must succeed. That probability is \(1-(0.60)^2=0.64\). Since \(P(K_1=2)=(0.40)^2=0.16\), this path contributes \(0.16(0.64)=0.1024\). Adding the disjoint paths gives the exact success probability \(0.0768+0.1024=0.1792\).
The trial enrolls two patients with probability 1 and enrolls two additional patients with probability \(P(K_1>0)=1-(0.60)^2=0.64\). Therefore its expected sample size is \(2+2(0.64)=3.28\) patients.
Suppose a simulation of 50,000 trials reports a success rate of 0.1790 and a mean sample size of 3.276. For the success rate, the benchmark Monte Carlo standard error is \(\sqrt{0.1792(0.8208)/50{,}000}=0.001715\); the difference, \(-0.0002\), is about \(-0.12\) standard errors. For sample size, each trial has \(N=2\) with probability 0.36 and \(N=4\) with probability 0.64. Thus \(\operatorname{Var}(N)=4(0.36)+16(0.64)-(3.28)^2=0.9216\), and the standard error of the simulated mean is \(\sqrt{0.9216/50{,}000}=0.004293\). The observed difference, \(3.276-3.28=-0.004\), is about \(-0.93\) standard errors. Both outputs are compatible with the exact results.
Use Limiting and Boundary Checks Carefully
Exact calculations are not always convenient, but simple limiting and boundary cases can still expose errors. Under a Bernoulli model with \(p=0\), every generated outcome should be zero; with \(p=1\), every outcome should be one. A trial rule with a futility threshold should behave deterministically when every possible interim outcome is known. These checks are especially useful for detecting indexing mistakes, incorrect inequality directions, and mishandled stopping conditions.
A limiting check should have a mathematical basis, not just an expectation that “more data looks stable.” For independent observations with equal finite variance, the variance of their mean decreases as \(\sigma^2/n\); with unequal variances, use the full sum-of-variances formula above. A decreasing variance does not mean every simulated sequence must move smoothly toward its limit. Random fluctuation remains at every finite sample size.
For a stochastic benchmark, a practical comparison reports the target, the simulated estimate, the number of replications, and the Monte Carlo standard error. A rule such as flagging a discrepancy larger than three Monte Carlo standard errors can be a useful screening convention, but it is not a universal pass/fail theorem. Under an approximately normal error distribution, a three-standard-error discrepancy is unusual for one prespecified check, yet checking many quantities makes at least one unusual result more likely. Investigate patterns and model assumptions rather than relying on a single threshold.
Diagnostics, Alternatives, and Limitations
When simulation and benchmark disagree, first verify that they refer to the same quantity, scale, and conditioning event. Then check the stated assumptions and calculation: Are outcomes independent where required? Does the benchmark use the same parameter values? Does the adaptive rule stop at the specified time? Are all paths included exactly once? For a rate, are there \(R\) independent replications, or did shared random numbers create dependence? A mislabeled denominator or an off-by-one stage boundary can create a discrepancy that looks like a statistical problem.
Three complementary approaches have different strengths:
| Approach | Useful for | Main limitation |
|---|---|---|
| Hand calculation or exact enumeration | Small discrete examples, exact probabilities, and rule paths | May become impractical as the state space grows |
| Analytical moments or known distributions | Checking generators, estimators, and summaries under tractable models | Only checks the features covered by the result and its assumptions |
| Limiting and boundary cases | Finding direction, threshold, and edge-case errors | Correct behavior at a few extremes does not establish correctness elsewhere |
No benchmark covers every combination of model, parameter, and adaptive path. A component can match its mean and variance while still having the wrong distribution, and an implementation can pass a boundary check while failing at an interior value. Use multiple benchmarks aimed at distinct failure modes. These checks support confidence in the implementation; they do not replace specification review or establish that the model is appropriate for the clinical question.
Common Mistake
Another frequent error is to treat a Monte Carlo standard error as a tolerance that guarantees correctness. It quantifies random error under the assumed simulation model; it cannot excuse a mismatch in the model, dependence between replications, or a wrong decision rule. A discrepancy within one standard error does not certify the code, just as a moderately large discrepancy does not by itself identify its cause.
Specify the probability law, parameter values, independence assumptions, and precise output being checked.
Use an exact probability, a moment calculation, path enumeration, or a justified boundary or limiting result.
Report the replication count and an appropriate Monte Carlo standard error for the simulated summary.
Check assumptions, units, branch logic, denominators, and whether the mismatch is larger or more systematic than Monte Carlo variation can reasonably explain.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Two independent observations have variances 3 and 7. What is the variance of their average?
- A simulated event rate is based on 10,000 independent replications, with exact event probability 0.10. What is its Monte Carlo standard error?
- Why is checking a simulated mean against its exact expectation not enough to establish that the full distribution is correct?
- For independent observations with variances \(\sigma_i^2\), state the condition expressed by the variance formula that ensures the variance of their average tends to zero as \(n\) increases.
Solutions. (1) The average has variance \((3+7)/2^2=10/4=2.5\). (2) The standard error is \(\sqrt{0.10(0.90)/10{,}000}=\sqrt{0.000009}=0.003000\). (3) Different distributions can share the same mean, so additional checks such as variance, event probabilities, or distributional features are needed. (4) The required condition is \(\sum_{i=1}^n\sigma_i^2/n^2\to0\); independence alone does not imply it, while uniformly bounded variances are sufficient.