Monte Carlo Error in Resampling Inference

Why Resampling Results Have Monte Carlo Error

A bootstrap distribution or permutation distribution is usually not calculated in full. Instead, an algorithm generates a finite number of resamples and uses their statistics to approximate a target quantity. Even with the observed data held fixed, a different sequence of random resamples would generally give a slightly different quantile or p-value. This variation is Monte Carlo error; it is distinct from the sampling uncertainty that the resampling procedure is intended to assess.

Let $G$ denote the conditional distribution of a resampled statistic, given the observed data, and let $B$ be the number of independently generated resamples. The simulation error typically decreases at rate $B^{-1/2}$. The constant in that rate matters: estimating an extreme quantile or a small tail probability can require many more replicates than estimating a central quantity. A reported result should therefore be accompanied by enough information to judge its numerical precision.

Definition: For a resampling target $\theta$, the Monte Carlo error is the randomness in its simulation-based estimate conditional on the observed data and the resampling scheme. Its Monte Carlo standard error (MCSE) is the standard deviation of that estimate over repeated runs of the simulation, with the data and scheme fixed.

Monte Carlo Error for a Bootstrap Quantile

Suppose $T_1^*,\ldots,T_B^*$ are independent draws from a bootstrap distribution with cumulative distribution function $G$. For $0<p<1$, define its $p$-quantile as $q_p=\inf\{x:G(x)\ge p\}$. The simulated quantile $\hat q_p$ is the corresponding quantile of the $B$ generated values. Under continuity and positive density at the target, the quantile's MCSE has a simple approximation.

Theorem: Suppose $G(q_p)=p$, and $G$ has a continuous density $g$ in a neighborhood of $q_p$, with $g(q_p)>0$. If $T_1^*,\ldots,T_B^*$ are independent draws from $G$, then $$ \sqrt{B}(\hat q_p-q_p)\xrightarrow{d}N\left(0,\frac{p(1-p)}{g(q_p)^2}\right). $$ Thus the approximate Monte Carlo standard error is $\sqrt{p(1-p)}/(\sqrt{B}\,g(q_p))$.

To justify the formula, let $\hat G_B$ be the empirical distribution function of the resamples. The empirical quantile satisfies $\hat G_B(\hat q_p)$ approximately equal to $p$. A first-order expansion around $q_p$, using differentiability of $G$, gives the quantile representation $$ \hat q_p-q_p=-\frac{\hat G_B(q_p)-p}{g(q_p)}+o_p(B^{-1/2}). $$ The empirical-CDF fluctuation at $q_p$ is the average of $B$ independent Bernoulli variables with success probability $p$. The central limit theorem therefore gives $$ \sqrt{B}\{\hat G_B(q_p)-p\}\xrightarrow{d}N(0,p(1-p)). $$ Substitution into the representation proves the stated limit and its variance. In practice, $g(q_p)$ is unknown and must be estimated, for example by a local density estimate from the resampled statistics. A poor density estimate can make this MCSE unreliable.

Worked Example: Monte Carlo Error for a Bootstrap Quantile

Suppose a bootstrap statistic has a distribution close to $N(0,1)$, and the goal is its 95th percentile. For a standard normal distribution, $q_{0.95}\approx1.64485$, and its density there is approximately $g(q_{0.95})=0.1031$. With $B=11{,}200$ bootstrap draws, the estimated MCSE is $$ \frac{\sqrt{0.95(0.05)}}{\sqrt{11{,}200}(0.1031)} = \frac{0.21794}{(105.83)(0.1031)} \approx 0.01997. $$

Thus a reported simulated 95th percentile has Monte Carlo variation of about two hundredths on the statistic's scale, under the smooth-distribution approximation. To target an MCSE no greater than $0.02$, solve the formula for $B$: $$ B\ge \frac{p(1-p)}{g(q_p)^2(0.02)^2} =\frac{0.0475}{(0.1031)^2(0.0004)} \approx 11{,}172. $$ Rounding up and using $11{,}200$ draws meets this approximate target. This calculation concerns Monte Carlo precision, not the sampling precision of the bootstrap confidence limit itself.

The density-based formula requires a smooth distribution with positive density near the quantile. Resampling distributions can instead be discrete, have ties, or be too irregular for a stable local density estimate. In those cases, order-statistic intervals give a direct alternative. For continuous $G$, the number of simulated values below the true $p$-quantile has a $\operatorname{Binomial}(B,p)$ distribution. Consequently, suitably chosen order statistics of $T_1^*,\ldots,T_B^*$ bracket $q_p$ with a specified probability. This interval describes Monte Carlo uncertainty in the quantile; its width can be reduced by increasing $B$.

Monte Carlo Error for a Permutation p-Value

For a fixed data set and a specified permutation scheme, let $\pi$ be the probability that a randomly drawn allowed permutation produces a statistic at least as extreme as the observed statistic. This is the tail probability under the permutation distribution. If $B$ permutations are drawn independently with replacement and $K$ are at least as extreme, then $K$ has a $\operatorname{Binomial}(B,\pi)$ distribution. The direct estimate $K/B$ has conditional variance $\pi(1-\pi)/B$, so its MCSE is estimated by replacing $\pi$ with $K/B$.

Result: For $K\sim\operatorname{Binomial}(B,\pi)$, the estimator $\hat\pi=K/B$ is unbiased and $$ \operatorname{Var}(\hat\pi)=\frac{\pi(1-\pi)}{B}, \qquad \widehat{\operatorname{MCSE}}(\hat\pi) =\sqrt{\frac{\hat\pi(1-\hat\pi)}{B}}. $$ The variance follows because $K$ is a sum of $B$ independent Bernoulli indicators, each with variance $\pi(1-\pi)$; dividing their sum by $B$ divides its variance by $B^2$.

A common Monte Carlo test reports $(K+1)/(B+1)$, rather than $K/B$. The added one corresponds to including the observed arrangement alongside the $B$ sampled permutations. This avoids a reported p-value of zero and, when the observed arrangement and sampled arrangements are exchangeable under the null, gives a valid finite-simulation test; ties handled by counting statistics at least as extreme make the test conservative. The plus-one value is a test p-value, not the unbiased estimate of the fixed tail probability $\pi$. Its simulation uncertainty is still of binomial order $B^{-1/2}$; for the plus-one estimator its conditional variance is $B\pi(1-\pi)/(B+1)^2$.

Worked Example: Precision of a Simulated Permutation p-Value

A test samples $B=9{,}999$ permutations and finds $K=37$ statistics at least as extreme as the observed statistic. The add-one Monte Carlo p-value is $$ \hat p_{\mathrm{MC}}=\frac{K+1}{B+1} =\frac{38}{10{,}000}=0.0038. $$

The direct estimate of the permutation tail probability is $37/9{,}999\approx0.00370037$. Its estimated MCSE is $$ \sqrt{\frac{(0.00370037)(1-0.00370037)}{9{,}999}} \approx 0.000607. $$ The uncertainty is substantial relative to the estimated tail probability: one should not interpret $0.0038$ as an exact numerical value of the permutation tail. More permutations can sharpen the simulation estimate. If the purpose is a test decision, the add-one p-value and the simulation design should be reported; if the purpose is to estimate $\pi$, report the tail estimate and a binomial interval as well.

Choosing a Replication Count

For a permutation tail probability, a normal approximation suggests that to obtain an MCSE no greater than $\varepsilon$, choose $$ B\ge \frac{\pi(1-\pi)}{\varepsilon^2}. $$ For a desired confidence margin $m$ at approximate two-sided confidence level $1-\alpha$, replace $\varepsilon$ by $m/z_{1-\alpha/2}$. Since $\pi$ is usually unknown before simulation, a pilot run can provide a planning value. If no credible planning value is available, the bound $\pi(1-\pi)\le 1/4$ gives a conservative count.

Worked Example: Planning Permutation Replicates

Suppose a pilot suggests a permutation tail probability near $0.01$, and the desired approximate 95% margin of error is $0.002$. Using $z_{0.975}=1.96$, the normal-approximation planning calculation is $$ B\ge \frac{(1.96)^2(0.01)(0.99)}{(0.002)^2} =\frac{3.8416(0.0099)}{0.000004} =9{,}507.96. $$ Thus at least $9{,}508$ draws are indicated by this approximation. If the pilot estimate is not trustworthy, using the worst-case variance $1/4$ instead gives $$ B\ge\frac{(1.96)^2(0.25)}{(0.002)^2}=240{,}100. $$

These are planning approximations, not guarantees of interval coverage. Near zero or one, a normal approximation to a binomial proportion may be poor; an exact binomial interval or a suitable score interval is preferable for assessing uncertainty. The same principle applies to bootstrap quantiles: estimate the MCSE under stated assumptions, then increase $B$ until it is small relative to the precision relevant to the analysis.

Stopping After Repeated Inspection

A fixed replication count chosen in advance allows ordinary fixed-sample confidence intervals. A different issue arises when the analyst repeatedly checks the estimate and stops as soon as its ordinary confidence interval looks sufficiently narrow. A 95% interval computed at each check does not generally retain 95% coverage at the random stopping time: repeated opportunities to stop when an interval happens to miss the target increase the overall error probability.

For a permutation tail probability, a simple remedy is an alpha-spending bound. Let $I_i$ indicate whether the $i$th independently drawn permutation is at least as extreme as the observed statistic, and let $\hat\pi_t=t^{-1}\sum_{i=1}^t I_i$. Choose a total error probability $\alpha$ and set $$ \alpha_t=\frac{6\alpha}{\pi^2t^2}, \qquad r_t=\sqrt{\frac{\log(2/\alpha_t)}{2t}}. $$ Hoeffding's inequality gives $P(|\hat\pi_t-\pi|>r_t)\le\alpha_t$. Since $\sum_{t=1}^{\infty}\alpha_t=\alpha$, the union bound proves $$ P\left(\text{there exists }t\ge1: |\hat\pi_t-\pi|>r_t\right) \le\sum_{t=1}^{\infty}\alpha_t =\alpha. $$ Therefore the intervals $[\hat\pi_t-r_t,\hat\pi_t+r_t]$, intersected with $[0,1]$, cover $\pi$ simultaneously at every inspection time with probability at least $1-\alpha$. Stopping when their width is below a preselected tolerance preserves this guarantee. The bound is conservative, but its validity does not depend on when the computation stops.

Worked Example: A Time-Uniform Stopping Bound

Take $\alpha=0.05$ and inspect the permutation estimate after $t=140{,}000$ independent draws. The allocated error probability is $$ \alpha_t=\frac{6(0.05)}{\pi^2(140{,}000)^2} \approx 1.55\times10^{-12}. $$ Then $$ \log(2/\alpha_t)\approx 27.89, \qquad r_t=\sqrt{\frac{\log(2/\alpha_t)}{280{,}000}} \approx\sqrt{\frac{27.89}{280{,}000}} \approx0.00998. $$

At this inspection time, the time-uniform interval has radius about $0.00998$. Its conservatism is visible: even after many draws, it may be wider than the desired precision for a small p-value. The guarantee is nevertheless valid if the analyst has inspected the estimate repeatedly and stops based on the interval. A sharper confidence sequence may improve efficiency, but it must itself provide time-uniform coverage; repeatedly applying a fixed-sample interval is not a substitute.

The same general idea applies to a bootstrap quantile, but the sequential target is not a Bernoulli mean. One can construct a simultaneous confidence band for the bootstrap CDF $G$, for example using the Dvoretzky–Kiefer–Wolfowitz inequality with an alpha-spending sequence, and invert that band to obtain a confidence set for $q_p$. Stopping when the inverted set is sufficiently narrow then retains the simultaneous coverage guarantee. Alternatively, use fixed $B$ and an order-statistic interval. A density-based quantile MCSE can guide planning, but repeatedly checking a pointwise MCSE alone does not establish sequential coverage.

Common Pitfalls and Practical Reporting

Confusing two kinds of uncertainty. Monte Carlo error describes finite simulation conditional on the data; it does not replace the sampling uncertainty represented by the bootstrap or permutation analysis.
Reporting a small p-value with false precision. A result such as $0.0038$ based on only a few dozen exceedances has appreciable simulation uncertainty. Report the number of draws and exceedances, and increase the simulation count if the distinction matters.
Using the quantile formula without its conditions. The density-based MCSE requires a positive, smooth density at the target quantile. For discrete or irregular bootstrap distributions, use order-statistic methods instead.
Stopping on an ordinary interval after repeated checks. Its nominal coverage applies at a fixed sample size, not automatically at a data-dependent stopping time. Use a prechosen replication count or a valid time-uniform procedure.

Key takeaway: Monte Carlo precision is a design choice. Quantile and tail-probability estimates typically improve at rate $B^{-1/2}$, but their constants differ. Set a precision target, use a method whose assumptions fit the resampling distribution, and account for repeated inspection whenever the simulation can stop adaptively.

Check Your Understanding

Try these questions before moving on.

Under what regularity conditions does the density-based MCSE formula for a bootstrap quantile apply?
Why can an order-statistic interval be preferable to estimating the density at a bootstrap quantile?
If $K$ of $B$ independently sampled permutations are at least as extreme as observed, what is the variance of $K/B$ conditional on the data?
For an estimated permutation tail probability near $0.01$, what replication-count formula gives an approximate target margin of error?
Why can repeatedly checking an ordinary 95% confidence interval invalidate its nominal coverage at the stopping time?
How does an alpha-spending sequence make repeated permutation checks compatible with an overall error guarantee?