Controlling Numerical Error in Small P-Values

Why Small Sequential P-Values Need Special Checks

A sequential p-value can be a small probability under the null hypothesis: the probability of an outcome at least as extreme as the observed stopped outcome, according to a specified ordering. Earlier material on stagewise and likelihood-ratio orderings defines which outcomes belong to that tail. Once the ordering and stopping rule are fixed, the numerical task is to evaluate the probability of that region under the null model.

The challenge is that a computation can look stable while still being inaccurate. A tiny omitted tail may be important relative to a tiny p-value. A quadrature routine may miss a narrow region or lose precision when adding very small quantities. Ordinary Monte Carlo may produce few or no tail events, leaving substantial uncertainty. These are distinct problems and require different checks.

Definition: Let \(Y\) denote the vector of statistics or increments used to represent a stopped sequential outcome, and let \(R_{\mathrm{obs}}\) be the region deemed at least as extreme as the observed outcome under the chosen ordering. The sequential p-value is \(p=P_{H_0}(Y\in R_{\mathrm{obs}})\). Numerical error is error in computing this probability for the specified null model and region; it is not a change in the ordering, stopping rule, or statistical model.

The calculations below use a simple two-stage illustration. Under \(H_0\), suppose the stagewise increments \(U_1\) and \(U_2\) are independent, each with distribution \(N(0,1)\). The event of interest is that the first increment is below zero and the second exceeds 5. This example is deliberately simple so the probability can also be checked analytically. In an actual sequential analysis, the joint distribution and tail region must match the design, the stopping rule, and the ordering being used.

Three Sources of Numerical Error

Truncation error occurs when an unbounded integral is evaluated only over a finite range. It can be controlled by increasing the range and, preferably, bounding the omitted probability analytically.

Quadrature error is error from approximating an integral with a finite numerical rule. Refining the grid, increasing the quadrature order, or tightening adaptive tolerances can reveal sensitivity, but agreement between refinements alone is not a proof of accuracy.

Monte Carlo error arises because simulation estimates a probability using a finite random sample. It persists even if the simulation code correctly generates the null distribution and checks the right event. For a rare event, a very large simulation may still yield only a small number of tail observations.

These computational errors should be distinguished from an incorrect covariance matrix, an incorrectly implemented stopping rule, or an event region that does not represent the chosen ordering. More digits do not correct a wrong model or event. A reliable check starts by verifying those inputs, then assesses the numerical method.

Worked Example: Bounding Error From a Truncated Tail

For the illustrative event \(U_1<0,\ U_2>5\), independence gives

p=P(U_1<0)P(U_2>5)=\frac{1}{2}\{1-\Phi(5)\}.

Using \(1-\Phi(5)\approx 2.8665157\times10^{-7}\), the probability is

p\approx \frac{1}{2}(2.8665157\times10^{-7}) =1.43325785\times10^{-7}.

Suppose numerical integration includes only values \(U_2\leq6\). The omitted probability is

P(U_1<0,\ U_2>6) =\frac{1}{2}\{1-\Phi(6)\} \approx\frac{1}{2}(9.865876\times10^{-10}) =4.932938\times10^{-10}.

Thus the truncated integral is about \(1.4283249\times10^{-7}\), below the full probability by \(4.932938\times10^{-10}\). Relative to the full probability, the omission is about \(0.3441\%\). It may look small on an absolute scale, but it is not negligible if a sub-tenth-of-a-percent relative accuracy is needed.

A useful bound follows from the standard normal density \(\phi\). For \(c>0\), since \(x/c\geq1\) for \(x\geq c\),

1-\Phi(c)=\int_c^\infty \phi(x)\,dx \leq \frac{1}{c}\int_c^\infty x\phi(x)\,dx =\frac{\phi(c)}{c}.

At \(c=6\), this gives \(1-\Phi(6)\leq \phi(6)/6\approx1.0127\times10^{-9}\). After multiplying by \(1/2\), the omitted part of this sequential event is at most \(5.063\times10^{-10}\). A conservative bound of this kind is preferable to assuming that a chosen integration limit is sufficiently far into the tail.

Quadrature: Refine, Rescale, and Cross-Check

For a sequential p-value written as a multivariate integral, the integrand may concentrate in a small part of the domain. A uniform grid can waste evaluations where the density is negligible and undersample the region that determines the tail probability. In addition to refining the integration rule, consider whether a change of variables can place the important tail region on a numerically convenient scale. If probabilities are extremely small, computing a log probability or using a tail-specific routine can avoid loss of precision from subtracting a number near one from one.

A practical convergence check compares results under meaningful refinements: for example, increasing quadrature order, tightening an adaptive tolerance, extending the integration range, and changing the parameterization. Record both the absolute difference and the relative difference. For two positive estimates \(p_1\) and \(p_2\), the relative difference using \(p_2\) as the reference is \(|p_1-p_2|/p_2\). The reference choice should be stated; a percentage without a denominator can be ambiguous.

Worked Example: Reading a Quadrature Refinement Check

For the same illustrative event, suppose an initial quadrature calculation reports \(1.4319\times10^{-7}\) and a refined calculation reports \(1.4331\times10^{-7}\). The analytic value above is \(1.43325785\times10^{-7}\), available here as an independent check.

The change between the two numerical calculations is

1.4331\times10^{-7}-1.4319\times10^{-7} =1.2\times10^{-10}.

Relative to the refined estimate, this change is \(1.2\times10^{-10}/(1.4331\times10^{-7})\approx0.000837\), or about \(0.0837\%\). The refined estimate differs from the analytic value by

1.43325785\times10^{-7}-1.4331\times10^{-7} =1.5785\times10^{-11},

which is about \(0.0110\%\) of the analytic value. In this particular check, refinement moves the result toward the independent value. In a real calculation, successive agreement would be encouraging but would not establish that the event, null distribution, or integration limits were correctly specified. A separate check of those inputs and, where practical, an independent computational method remain important.

Monte Carlo Error for Rare Events

With \(N\) independent simulations under \(H_0\), let \(K\) be the number that fall in \(R_{\mathrm{obs}}\). The Monte Carlo estimate is \(\hat p=K/N\). Under this simulation model, \(K\) has a binomial distribution with parameters \(N\) and \(p\), so the estimated standard error is

SE(\hat p)=\sqrt{\frac{\hat p(1-\hat p)}{N}}.

When \(p\) is small, the relative standard error is approximately \(1/\sqrt{Np}\). This approximation explains why a simulation that seems large in absolute terms may be poor for a very small tail probability: the expected number of tail events, \(Np\), may still be low. If no events are observed, the estimate is zero, but that does not establish that the true probability is zero.

Worked Example: Assessing a Rare-Event Simulation

Suppose \(N=100{,}000{,}000\) independent null simulations produce \(K=16\) outcomes in the tail region. Then

\hat p=\frac{16}{100{,}000{,}000}=1.6\times10^{-7}.

The estimated standard error is approximately

SE(\hat p) =\sqrt{\frac{(1.6\times10^{-7})(1-1.6\times10^{-7})}{10^8}} \approx4.0\times10^{-8}.

The relative standard error is about \(4.0\times10^{-8}/(1.6\times10^{-7})=25\%\). A rough normal-approximation interval is \(\hat p\pm1.96SE\), or about \([8.16\times10^{-8},\,2.384\times10^{-7}]\). Because this interval is based on only 16 events, a binomial or Poisson interval is preferable for a formal uncertainty report; the rough calculation already shows that the simulation does not support several precise digits.

For a target relative standard error near \(10\%\), the approximation \(1/\sqrt{Np}\approx0.10\) requires \(Np\approx100\). If \(p\) is about \(1.43\times10^{-7}\), this means \(N\) should be roughly \(100/(1.43\times10^{-7})\approx7.0\times10^8\) simulations. A \(5\%\) target requires about \(Np=400\), or approximately \(2.8\times10^9\) simulations at that probability. These are planning approximations, not guarantees about a particular run.

A Precision Check That Matches the Claim

The appropriate check depends on what is being reported. If the claim is only that \(p<10^{-5}\), an absolute error much smaller than \(10^{-5}\), or a conservative upper bound below \(10^{-5}\), may be sufficient. If the result is reported as \(p=1.4\times10^{-7}\), the computation needs enough relative accuracy to support the displayed digits. A statement such as \(p<10^{-6}\) is not justified by a simulation that observed zero events in a sample too small to rule out probabilities near \(10^{-6}\).

Use the following checks together where feasible:

Verify that the null distribution, covariance structure, stopping rule, and ordering all match the intended sequential p-value.
Increase integration limits and bound the omitted probability; do not rely only on a visual impression that the density is small at the boundary.
Refine quadrature settings and compare results on both absolute and relative scales. Check numerical stability, including under a different parameterization or implementation when possible.
For simulation, report \(N\), \(K\), the estimate, and an uncertainty measure. For very small probabilities, do not interpret zero observed events as a zero p-value.
Round only after the error assessment. Report no more digits than are supported by the numerical checks, and distinguish a computed estimate from a guaranteed bound.

These methods have complementary strengths. Quadrature can be efficient and reproducible when the integral is low-dimensional and well behaved, but it can struggle with high-dimensional regions or complicated boundaries. Direct Monte Carlo is flexible and easy to validate conceptually, but its relative error grows rapidly as the probability decreases. Tail bounds can certify that a discarded region is negligible, but only when the bound is valid for the model and event. Importance sampling can target rare regions more efficiently, but requires careful weighting and a separate assessment of its variance and implementation.

Common Pitfalls and Limitations

Common mistake: “The absolute difference between two calculations is tiny, so the small p-value is accurate.”

Correction: Compare the difference with the p-value itself as well as with the required reporting precision. A difference that is tiny in ordinary units can be a large fraction of a probability near zero. Also check the tail omitted by truncation and verify that both calculations evaluate the same event under the same null model.

A second pitfall is treating repeated agreement from one numerical routine as independent confirmation. Two runs with the same truncation, parameterization, and implementation can share the same bias. Another is reporting a simulated estimate with many digits when only a handful of events were observed. More computation can reduce Monte Carlo uncertainty, but it cannot fix an incorrectly specified tail region or an inappropriate null covariance structure.

No numerical check eliminates all uncertainty. A truncation bound controls only the omitted part it covers; quadrature refinement is diagnostic rather than a universal error guarantee; and simulation intervals quantify random error under the simulation model, not coding or model error. The final reported precision should reflect the weakest relevant check.

Check Your Understanding

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

If a probability estimate changes from \(4.0\times10^{-8}\) to \(4.4\times10^{-8}\), what is the change relative to the refined estimate?
Why does truncating an integral at a point where the density is small not by itself show that the omitted probability is negligible?
A simulation uses \(N=10^7\) draws and observes \(K=0\) tail events. What can and cannot be concluded?
Name two checks, beyond simply rerunning the same quadrature routine, that can strengthen an accuracy assessment.

Solutions. (1) The absolute change is \(0.4\times10^{-8}=4.0\times10^{-9}\). Relative to the refined estimate, it is \((4.0\times10^{-9})/(4.4\times10^{-8})\approx0.0909\), or about \(9.09\%\). (2) A small density at the boundary does not control the total mass beyond it; the remaining tail must be integrated or bounded. (3) The estimate is zero, but the true probability is not thereby shown to be zero. The run gives limited information about a rare event; an appropriate binomial upper confidence bound can quantify what probabilities remain compatible with zero observed events. (4) Examples include extending the integration domain and bounding the omitted tail, changing quadrature order or parameterization, comparing with an independent method, and verifying the event region and null distribution against the design.)