Monte Carlo Error In Prior Predictive Summaries

Why simulation error matters

A prior predictive check asks what trial outcomes are plausible under the specified prior and sampling model. Often those outcomes are simulated rather than calculated exactly. The resulting summaries—such as a tail probability or a high quantile—are estimates based on a finite number of simulated trial replicates. They can change if the simulation is run again, even when the prior, trial design, and code are unchanged.

That run-to-run variation is Monte Carlo error. It is not uncertainty about the treatment effect in a particular trial, nor does it represent uncertainty about whether the prior is scientifically appropriate. It is numerical sampling variation introduced by using a finite simulation to approximate a prior predictive quantity. A useful check therefore reports not only a summary but also whether the number of replicates gives that summary adequate simulation precision.

For the usual calculations below, each replicate must be generated independently from the full prior predictive distribution: draw its parameter values from the prior, then draw its trial outcome from the sampling model conditional on those values. The resulting replicate-level summaries must also be independent and identically distributed. Reusing the same parameter draw across multiple outcomes, or otherwise inducing dependence among replicates, invalidates the simple standard-error formulas unless that dependence is explicitly accounted for.

Monte Carlo error for a tail probability

Each indicator \(I(T_m\geq c)\) is a Bernoulli random variable with probability \(p\). Thus \(K\sim\operatorname{Binomial}(M,p)\), \(E(\hat{p}_{MC})=p\), and \(\operatorname{Var}(\hat{p}_{MC})=p(1-p)/M\). Replacing the unknown \(p\) by \(\hat{p}_{MC}\) gives the estimated standard error above. This calculation concerns the precision of the simulated estimate of the probability, not the variability of \(T\) across hypothetical trials.

A normal-approximation interval is \(\hat{p}_{MC}\pm1.96\,\widehat{SE}(\hat{p}_{MC})\), but it can perform poorly when the event is rare or the count \(K\) is small. For a reportable interval, a binomial interval such as the Wilson interval is generally preferable in those cases. Regardless of interval method, simulation precision should be judged against a tolerance chosen for the intended decision, rather than by whether a rounded number looks stable.

Quantiles for discrete and continuous summaries

A prior predictive summary may be a count, a binary outcome, or another discrete quantity. Its distribution can have jumps, so a quantile should not be defined by requiring that its cumulative distribution equal the probability level. Let \(F(t)=P(T\leq t)\) be the cumulative distribution function of \(T\). For \(0<\tau<1\), use the generalized inverse definition:

For example, if \(T\) is Bernoulli with \(P(T=1)=0.2\), then \(F(0)=0.8\). Its median under the generalized inverse definition is \(q_{0.5}=0\), even though \(F(0)\ne0.5\). This matters when auditing trial counts or other discrete summaries: the appropriate quantile may lie at a jump, and an equation \(F(q_\tau)=\tau\) may have no solution.

For simulated values, define the empirical distribution function \(\hat{F}_M(t)=M^{-1}\sum_{m=1}^M I(T_m\leq t)\). The corresponding estimated quantile is \(\hat{q}_\tau=\inf\{t:\hat{F}_M(t)\geq\tau\}\), equivalently the order statistic \(T_{(k)}\) for \(k=\lceil M\tau\rceil\), where the values are sorted from smallest to largest. Ties are retained. This estimate changes in steps as the simulation count changes, particularly for discrete outcomes.

Quantifying Monte Carlo error in a quantile

For a continuous distribution with a positive, smooth density \(f\) near \(q_\tau\), a large-sample approximation to the quantile estimator's Monte Carlo standard error is \[ SE(\hat{q}_\tau)\ \approx\ \sqrt{\frac{\tau(1-\tau)}{M f(q_\tau)^2}}. \] This approximation explains why extreme quantiles and regions with low density can require many more replicates than central quantiles. The density in the formula is usually unknown and must be estimated, so the calculated standard error is itself approximate.

A short justification follows from the empirical distribution at the true quantile. For a continuous distribution, \(F(q_\tau)=\tau\), and the number of simulated values at or below \(q_\tau\) has a binomial distribution with parameters \(M\) and \(\tau\). Consequently, the empirical CDF at \(q_\tau\) has variance \(\tau(1-\tau)/M\). Near \(q_\tau\), a small change in the CDF is approximately \(f(q_\tau)\) times the corresponding change in the quantile. Dividing the CDF-scale standard deviation by \(f(q_\tau)\) gives the stated approximation. The argument requires continuity and positive density near the quantile; it is not a general formula for quantiles at jumps.

When the summary is discrete or has a jump near the target quantile, do not apply this density formula mechanically. Instead, assess the quantile using the simulated empirical distribution and uncertainty for its CDF. At a proposed value \(t\), the count of simulated values at or below \(t\) is binomial with probability \(F(t)\). Binomial confidence bounds for \(F(t)\), evaluated across relevant values, can indicate how well the simulation locates the first value where the CDF reaches \(\tau\). For a discrete outcome, report the quantile and the adjacent CDF values when useful; they show whether the quantile is determined by a large jump or is close to a boundary.

Choosing a replicate count for a precision goal

Adequacy is a practical decision: specify how much simulation error is acceptable for the summary and the use to which it will be put. For a tail probability, a target approximate 95% half-width \(h\) suggests \[ 1.96\sqrt{\frac{p(1-p)}{M}}\leq h, \qquad M\geq\frac{1.96^2p(1-p)}{h^2}. \] Because \(p\) is not known before simulation, a pilot estimate can guide the calculation. If a conservative planning bound is needed, use \(p(1-p)\leq0.25\), giving \(M\geq1.96^2(0.25)/h^2\). This is conservative for probabilities far from 0.5. For quantiles, use the density-based approximation only when its assumptions are plausible; otherwise examine the empirical CDF and quantify uncertainty at candidate values.

Diagnostics, alternatives, and limitations

A stable-looking estimate is not by itself evidence of adequate simulation. Check the event count for tail probabilities, the number and pattern of ties for quantiles, and whether the simulated replicates were generated independently. For quantiles, compare results across larger independent batches or a longer run and inspect the empirical CDF near the target level. Such comparisons are diagnostics, not substitutes for a precision calculation. If a local density estimate changes substantially with its bandwidth or nearby range, the density-based quantile standard error is not dependable.

If an exact prior predictive calculation is available, it can avoid Monte Carlo error for the quantity calculated exactly. Deterministic numerical integration may also be useful when the parameter space and likelihood permit it. Simulation remains flexible for complex trial designs, but finite simulation error must be reported and controlled. These methods address numerical precision; none determines whether the prior itself is scientifically credible.

The binomial formulas also rely on independent replicate-level indicators. If draws are dependent, the usual standard error can understate or overstate the actual error. Likewise, a pilot estimate can be a poor guide for a very rare event if it observes few or no occurrences. For such cases, report the event count and use an interval that respects the binomial structure, or increase the simulation effort. Do not interpret an estimated probability of zero as proof that the event is impossible under the prior predictive model.

Common mistake