Numerical Integration for Predictive Probability

Why Numerical Integration Is Useful

A posterior predictive probability averages future success probabilities over uncertainty about the treatment effect. In a simple normal model, that average can sometimes be calculated exactly. In more complicated success rules, it generally requires numerical integration. Quadrature provides a deterministic way to approximate the integral, and the normal model gives useful exact results for checking whether the implementation is correct.

The method is especially useful when a success rule depends on more than one future statistic, or when the posterior distribution is a mixture. Simulation, covered in the earlier material on Monte Carlo estimation of predictive probability, is flexible but has sampling error. Quadrature avoids that particular source of error, although it can still be inaccurate if the integration grid or tail calculation is poorly chosen.

Model and Integral

Suppose \(\theta\) is the treatment effect, and the interim data \(D\) give a posterior distribution \(\theta\mid D\sim N(m,v)\), where \(m\) is the posterior mean and \(v\) is the posterior variance. Let \(Y\) be a future estimate of the effect. Assume that, conditional on \(\theta\) and the interim data, \(Y\sim N(\theta,\tau^2)\), with \(\tau^2\) known and the future data independent of the interim data given \(\theta\). For a one-sided success rule \(Y>c\), success means that the future estimate exceeds the specified boundary \(c\), in the same units as \(\theta\) and \(Y\).

For the stated model, the conditional probability is \(\Phi((\theta-c)/\tau)\). Thus the predictive probability is an integral over the posterior. Standardizing the posterior with \(Z=(\theta-m)/\sqrt v\), where \(Z\sim N(0,1)\), gives an integral with respect to the standard normal density \(\phi\):

This standardization is convenient because the integration variable and weights do not depend on the posterior mean or variance. For a general success event, replace the conditional probability inside the integral by the probability of that event given \(\theta\) and the interim data. The same posterior averaging principle applies.

Gauss–Hermite Quadrature

An \(r\)-point Gauss–Hermite rule approximates an integral weighted by a Gaussian density using a finite weighted sum. Let \(x_i\) and \(w_i\), for \(i=1,\ldots,r\), be nodes and weights for the convention \(\int_{-\infty}^{\infty}e^{-x^2}f(x)\,dx\approx\sum_{i=1}^r w_i f(x_i)\). Since a standard normal variable has the same distribution as \(\sqrt{2}X\) when \(X\) has density proportional to \(e^{-x^2}\), the predictive probability is approximated by:

The factor \(1/\sqrt{\pi}\) and the \(\sqrt{2v}\) in the function argument are both important. Using nodes or weights defined for a different quadrature convention without changing the formula can produce a plausible-looking but incorrect answer. A practical implementation evaluates the conditional success probability at each transformed node and takes the weighted average.

Worked Example: A Three-Point Quadrature Check

An Exact Benchmark and Its Justification

For the single normal threshold rule, the integral has a closed form. One way to see this is to write \(Y=\theta+\varepsilon\), where, conditional on \(\theta\) and \(D\), \(\varepsilon\sim N(0,\tau^2)\) and is independent of \(\theta\mid D\). Since \(\theta\mid D\sim N(m,v)\), the sum is normally distributed with mean \(m\) and variance \(v+\tau^2\). Consequently:

This identity is both a result and a diagnostic. For a single linear threshold under the stated normal assumptions, use the closed form rather than quadrature if only the probability is needed. Quadrature remains valuable as a general computational approach for more complex success events, and this simpler case lets you test code, node conventions, and posterior parameterization before applying them to a harder integral.

Worked Example: An Extreme Lower-Tail Probability

Mixture Posteriors and Componentwise Integration

If the posterior is a mixture, integrate over each component and then average the component-specific probabilities using their posterior weights. Suppose \(\theta\mid D\) has components \(N(m_j,v_j)\) with weights \(\pi_j\), where \(\pi_j\geq0\) and \(\sum_j\pi_j=1\). For the normal future estimate and threshold above, each component has its own normal predictive distribution. The mixture predictive probability is:

Accuracy Checks and Failure Modes

A quadrature answer is an approximation, not a guarantee of accuracy. Check the calculation in ways that target different potential errors:

• Check parameterization. Confirm that \(v\) is the posterior variance and \(\tau^2\) is the conditional future variance. The normal integral uses \(\sqrt v\) to transform the posterior and \(\tau\) in the conditional probability.
• Check weights and nodes. For a normalized standard-normal quadrature rule, the weights used directly with \(\phi\) must sum to one. For the physicists’ Gauss–Hermite convention used above, the weights are instead divided by \(\sqrt{\pi}\).
• Increase the order. Compare results using progressively more quadrature points. Stability across orders is reassuring, although it does not by itself prove accuracy if the integrand is poorly resolved.
• Use an exact special case. In the single-threshold normal model, compare with the closed form derived above. A discrepancy can reveal a coding error, an incorrect tail direction, or an inconsistent variance.
• Inspect the tails. For probabilities very near zero or one, tiny absolute errors can be large relative errors. Direct computation of a small probability as \(1-\Phi(a)\) for large positive \(a\) can lose precision through subtraction. Use a numerically stable survival-function or log-CDF calculation where available.

Fixed-node quadrature can perform poorly when the region contributing most of the probability is narrow or far into a tail and the chosen nodes do not resolve it. An adaptive integrator can allocate effort to that region, but it too needs appropriate tolerances and stable evaluation of the integrand. Truncating a standard-normal integral to \([-K,K]\) also omits posterior mass \(P(|Z|>K)\); a small omitted mass is not automatically negligible relative to an extremely small predictive probability. For such cases, check convergence with different transformations or integration strategies and report the numerical tolerance or stability checks used.

Quadrature, Simulation, and the Closed Form

The methods are complementary. Use the closed form when its assumptions and event definition apply. For a low-dimensional integral, quadrature often gives a precise answer with no Monte Carlo sampling variability. Simulation is often easier to extend to complex or high-dimensional predictive models, but its uncertainty should be quantified, especially for rare events. Whenever possible, compare quadrature and simulation in a tractable case before relying on either for a more complicated calculation.

Common Mistake

A frequent error is to treat the posterior mean as if it were the true treatment effect, calculate a conditional success probability at that one value, and label the result predictive. That calculation ignores posterior uncertainty. The predictive probability instead averages conditional probabilities over the full posterior distribution; in the normal threshold example, this adds \(v\) to the future variance. A second error is to assume that a quadrature result is accurate simply because it is a deterministic number. Determinism removes simulation noise, not numerical integration error. Check the result against an exact benchmark where available and assess convergence, particularly in extreme tails.