Analytic Benchmarks For Simulation Verification

Why use analytic benchmarks?

A trial simulation can produce plausible-looking results even when its data generator, decision rule, or summary code contains an error. In simple settings, probability theory supplies results that can be calculated independently of the simulation. Comparing simulated outputs with those results is a practical way to check whether the simulation behaves as expected before using it for more complicated trial scenarios.

An analytic benchmark is not a requirement that a finite simulation match a theoretical value exactly. Simulated estimates vary by chance. The useful question is whether the difference is a reasonable amount of Monte Carlo variation, given the number of independent replicates. A large or persistent discrepancy can point to a coding or implementation problem; a small discrepancy is not, by itself, evidence of one.

Definition: An analytic benchmark is a known probability, distribution, expectation, or other operating characteristic under a specified data-generating model and analysis rule. A simulation verification check compares a simulated estimate of that quantity with its analytic benchmark, allowing for the simulation's Monte Carlo uncertainty.

The benchmark is meaningful only if its assumptions match the simulated setting. For example, an exact binomial probability assumes independent Bernoulli outcomes with a common event probability. A normal-theory benchmark for a sample mean assumes independent normal observations with the stated mean and variance. If the simulation intentionally adds clustering, unequal probabilities, missingness, or another feature, the simple benchmark no longer describes that simulation.

Quantifying the expected simulation variation

Suppose \(R\) independent trial replicates each produce an indicator \(I_r\), equal to 1 when a specified event occurs and 0 otherwise. Let the event probability under the model be \(q\). Then the number of event replicates \(K=\sum_{r=1}^R I_r\) has a binomial distribution, and the simulated event rate \(\hat{q}=K/R\) has mean \(q\) and variance \(q(1-q)/R\).

K \sim \operatorname{Binomial}(R,q), \qquad E(\hat{q})=q, \qquad SE(\hat{q})=\sqrt{\frac{q(1-q)}{R}}.

This standard error describes the typical Monte Carlo variation in the estimated event rate, not variation between patients within a trial. The benchmark value \(q\) is used in the calculation because it is known for this verification exercise. In an application where \(q\) is unknown, the estimated rate may instead be used to estimate its Monte Carlo standard error.

A convenient comparison is the standardized difference \((\hat{q}-q)/SE(\hat{q})\). Under the binomial model, this difference is approximately standard normal when the expected numbers of event and non-event replicates are sufficiently large. It is a diagnostic scale, not a rule that a particular cutoff proves the code correct or incorrect. If the expected count is small, direct binomial probabilities or an exact binomial interval are preferable to relying on a normal approximation.

Worked benchmarks

Worked Example: Checking a Bernoulli Data Generator

Suppose a generator is intended to produce independent binary outcomes with event probability \(p=0.30\). It generates \(R=5{,}000\) independent datasets, each containing \(n=20\) outcomes. Across all datasets, the simulated event proportion is 0.302. Check the pooled proportion against its benchmark.

Under the stated assumptions, there are \(5{,}000 \times 20=100{,}000\) independent Bernoulli outcomes. The number of events therefore has a binomial distribution with mean \(100{,}000(0.30)=30{,}000\), and the pooled proportion has standard error \(\sqrt{0.30(0.70)/100{,}000}=0.001449\).

SE(\hat{p})=\sqrt{\frac{p(1-p)}{Rn}} =\sqrt{\frac{0.30(0.70)}{100{,}000}} =0.001449, \qquad z=\frac{0.302-0.300}{0.001449}=1.380.

The observed proportion is 0.002, or about 1.38 standard errors, above the target. That difference is compatible with ordinary simulation variation; it is not a reason to demand that the generator return exactly 0.300. This pooled check assumes all generated outcomes are independent and use the same \(p\). If outcomes are correlated within datasets, the pooled standard error above is not appropriate, even if each individual outcome has marginal event probability 0.30.

Worked Example: Checking a Binomial Decision Probability

In a deliberately simple trial scenario, each of 10 participants has an independent binary outcome with event probability \(p=0.20\). The analysis declares a positive decision when at least 4 participants have the event. Find the exact probability of a positive decision, then assess a simulation that reports a positive-decision rate of 0.1225 across \(R=20{,}000\) independent replicates.

Let \(X\) be the event count in one replicate. The assumptions give \(X\sim\operatorname{Binomial}(10,0.20)\). The decision probability is the upper tail:

q=P(X\geq 4) =1-\sum_{k=0}^{3}{10\choose k}(0.20)^k(0.80)^{10-k} =0.1209.

For clarity, the probabilities for counts 0 through 3 are approximately 0.1074, 0.2684, 0.3020, and 0.2013. Their sum is 0.8791, giving an upper-tail probability of \(1-0.8791=0.1209\). Across 20,000 independent replicates, the Monte Carlo standard error at this benchmark is \(\sqrt{0.1209(1-0.1209)/20{,}000}=0.002305\).

z=\frac{0.1225-0.1209}{0.002305}=0.694.

The simulated rate is about 0.69 standard errors above the exact probability, a modest difference consistent with Monte Carlo variation. This check can reveal, for example, an incorrectly implemented threshold such as \(X>4\) instead of \(X\geq4\), or a generator using the wrong event probability. It does not establish that a more complex trial model is correct; it checks this simple scenario and this specific decision calculation.

Worked Example: Checking a Normal-Theory Decision Probability

Consider a simple continuous-outcome scenario in which the observations within each replicate are independent and follow \(X_i\sim N(5,4)\). Thus the population mean is 5 outcome units and the population variance is 4 squared outcome units. Each replicate contains \(n=25\) observations, and the rule declares a positive decision if \(\bar{x}>5.8\). Find the analytic probability of a positive decision and compare it with a simulated rate of 0.0242 from 10,000 independent replicates.

For independent normal observations, the sample mean is exactly normal with mean \(\mu\) and variance \(\sigma^2/n\). Here its standard deviation is \(\sqrt{4/25}=0.4\) outcome units. Standardizing the decision threshold gives \(z=(5.8-5)/0.4=2\), so the exact positive-decision probability is the upper tail above 2.

\bar{X}\sim N\left(5,\frac{4}{25}\right), \qquad q=P(\bar{X}>5.8)=1-\Phi(2)=0.02275.

The Monte Carlo standard error of the simulated rate is \(\sqrt{0.02275(1-0.02275)/10{,}000}=0.001491\). Its standardized difference from the benchmark is \((0.0242-0.02275)/0.001491=0.973\).

SE(\hat{q})=\sqrt{\frac{0.02275(1-0.02275)}{10{,}000}} =0.001491, \qquad z=\frac{0.0242-0.02275}{0.001491}=0.973.

The simulated rate is within one Monte Carlo standard error of the normal-theory value. The check depends on both the data-generating and analysis assumptions: observations must be independent normal draws with variance 4, the analysis must compute the mean on the original outcome scale, and the decision must use the strict threshold \(\bar{x}>5.8\). A change in variance, sample size, or inequality changes the benchmark.

Why the normal benchmark has this form

The normal-theory result follows directly from the model. If \(X_1,\ldots,X_n\) are independent normal random variables with mean \(\mu\) and variance \(\sigma^2\), the sum of independent normal variables is normal. Scaling that sum by \(1/n\) gives a normal sample mean. Its mean is \(\mu\), and independence makes its variance the sum of the individual variances divided by \(n^2\):

\bar{X}=\frac{1}{n}\sum_{i=1}^n X_i,\qquad E(\bar{X})=\frac{1}{n}\sum_{i=1}^n\mu=\mu,\qquad \operatorname{Var}(\bar{X})=\frac{1}{n^2}\sum_{i=1}^n\sigma^2=\frac{\sigma^2}{n}.

Therefore, for a rule that declares an event when \(\bar{X}>c\), with \(c\) a fixed threshold, the analytic operating characteristic is \(P(\bar{X}>c)=1-\Phi((c-\mu)/(\sigma/\sqrt{n}))\). This result is exact under the independent normal model. For nonnormal observations, the same expression may be an approximation under additional conditions, but it is not an exact benchmark merely because a sample mean is being calculated.

A disciplined comparison

A useful verification comparison has three parts: state the probability model, derive the target quantity independently, and calculate how much a finite simulation could vary around that target. Then inspect discrepancies in the context of the code path being checked. For a single event rate, a difference near one or two Monte Carlo standard errors is unsurprising. A difference of many standard errors merits investigation, especially if it recurs under independent runs or grows with the number of replicates.

Do not treat a standardized difference as a mechanical pass/fail test. If many benchmarks are checked, some will appear unusually discrepant by chance. Also, a result close to the benchmark does not prove all code is correct: distinct errors can cancel, and a benchmark exercises only the model and outputs it covers. Rechecking with a second parameter setting, another summary, or an independently derived probability can make the evidence more informative.

Approach	What it provides	Main limitation
Analytic benchmark	Known value for a tractable model, with a calculable simulation error scale	Applies only when the benchmark assumptions match the simulation
Simulation without a benchmark	Estimates operating characteristics for models that may be too complex for exact calculation	Agreement with a plausible-looking result is hard to judge independently
Independent implementation	A second code path can expose implementation-specific mistakes	Two implementations may share the same mistaken assumption or formula

Diagnostics, failure modes, and limits

Check the unit of replication. If the operating characteristic is a trial-level decision rate, the independent units for its Monte Carlo standard error are independent trial replicates, not the participants pooled across trials.
Check the probability model. Clustering, unequal event probabilities, dependence across replicates, or reuse of random draws can invalidate the simple binomial standard error.
Check the rule’s boundary. “At least 4” means \(X\geq4\), whereas “more than 4” means \(X>4\). These are different probabilities for discrete data.
Check the scale and variance. A normal benchmark uses the variance in the model and the sample size in the analysis. Confusing variance with standard deviation changes the probability.
Check numerical and coding details. A tail probability may be lost through subtracting nearly equal rounded numbers; use a reliable probability calculation when tails are very small. Also confirm that simulated data are not silently rounded or transformed before the decision is applied.

Analytic benchmarks are deliberately narrow. Real trial simulations may use hierarchical models, adaptive rules, missing data, complex randomization, and posterior calculations for which no simple exact operating characteristic is available. In those settings, the benchmark examples still help test simplified configurations and individual outputs, but they cannot certify the full simulation. The appropriate benchmark must match the component being checked; otherwise, a disagreement can reflect a mismatched model rather than an implementation error.

Key takeaway: Compare simulated results with a known value and its Monte Carlo uncertainty, not with exact equality. A benchmark is informative only when its model, independence assumptions, analysis, and decision rule correspond to the simulated setting.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A simulation estimates an event probability of 0.40 from 1,000 independent replicates. What is the Monte Carlo standard error at the benchmark value 0.40?
For \(X\sim\operatorname{Binomial}(5,0.10)\), write the exact expression for \(P(X\geq2)\). Why is it not enough to compare a simulated rate with 0.10?
Independent observations have distribution \(N(10,9)\), and each replicate contains 36 observations. What is the distribution and standard deviation of the sample mean?
A simulated positive-decision rate differs from its analytic benchmark by 0.001, and its Monte Carlo standard error is 0.004. What can be concluded from that comparison?

Solutions. (1) \(SE=\sqrt{0.40(0.60)/1{,}000}=0.01549\). (2) \(P(X\geq2)=1-P(X=0)-P(X=1)=1-(0.90)^5-5(0.10)(0.90)^4\), approximately 0.08146. The decision event is a count threshold, so its probability is a binomial tail, not the individual event probability. (3) The mean is 10 and the variance is \(9/36=0.25\), so \(\bar{X}\sim N(10,0.25)\) with standard deviation 0.5 outcome units. (4) The discrepancy is \(0.001/0.004=0.25\) Monte Carlo standard errors, which is small and consistent with simulation variation; it does not prove the code is correct.