Why Numerical Stability Matters
Normalized prediction distribution errors (NPDEs) use simulated data to compare observations with a model’s predictive distribution. For a subject with repeated observations, the calculation also accounts for dependence among observations by estimating their covariance and transforming the resulting errors. This makes NPDEs useful for evaluating predictive performance, but it also creates numerical steps that can fail even when the model fit itself has completed successfully.
The key distinction is between a statistical problem and a computational problem. An observation that is unusual under the predictive distribution may be a meaningful model diagnostic. An infinite score caused by a simulated rank at the boundary, or a failed matrix factorization caused by redundant observations, is not automatically evidence of poor predictive performance. Diagnose the calculation before interpreting its output.
The marginal scores are intended to be approximately standard normal under the fitted predictive model. The covariance step addresses their dependence within a subject. With a correct predictive model and adequate simulation, the decorrelated scores should be approximately standard normal and uncorrelated. They are not guaranteed to be exactly independent or standard normal in a finite simulation. The earlier material on derivation, decorrelation, and finite-sample behavior of NPDEs supplies the statistical context; this tutorial focuses on computational reliability.
Finite Simulation and Extreme Ranks
A predictive percentile is often estimated by ranking an observed value among \(M\) simulated values. The empirical cumulative distribution function can attain 0 or 1 when the observation lies beyond every simulated value. Applying the standard normal quantile function \(\Phi^{-1}\) to either endpoint gives negative or positive infinity. This is a boundary problem in the empirical estimate, not evidence that the underlying predictive probability is exactly zero or one.
One common finite-sample convention uses a mid-rank plotting position. If \(r\) is the number of simulated values strictly below the observation and \(t\) is the number tied with it, define \( \hat F=(r+t/2+0.5)/(M+1) \). This convention keeps the estimate strictly between 0 and 1, including when the observed value is outside the simulated range. Other implementations use different plotting positions or smooth the empirical distribution; the convention should be documented and applied consistently.
Worked Example: Converting a Boundary Rank to a Finite Score
Suppose 99 simulated concentrations are generated at a sampling time and the observed concentration exceeds all 99. For this illustration, use the mid-rank plotting position above, with \(M=99\), \(r=99\), and \(t=0\).
The estimated percentile is \( \hat F=(99+0+0.5)/(99+1)=99.5/100=0.995 \). The corresponding normal score is \( \Phi^{-1}(0.995)\approx 2.576 \).
The score is large and positive, consistent with an observation above all simulated values, but finite. By contrast, the unsmoothed empirical CDF gives \(99/99=1\) if defined as the proportion less than or equal to the observation, and \( \Phi^{-1}(1)=+\infty \). The finite plotting position prevents an infinite computer value; it does not make the observation less extreme or add information beyond the 99 simulations.
Because the observation lies beyond every simulated value, the estimate is especially sensitive to \(M\) and to the plotting-position rule. Report the simulation count and rule, and consider increasing \(M\) when boundary ranks occur frequently or materially affect the diagnostic.
Ties need attention as well. They can arise from rounded measurements, censoring conventions, discrete outcomes, or simulations with repeated values. An implementation that ignores ties may assign arbitrary ranks to identical values. A documented mid-rank rule, or a suitable randomized tie procedure when justified, avoids treating tied values as if their ordering were meaningful. The rule must match the outcome and predictive distribution; adding arbitrary jitter simply to eliminate ties can distort the check.
Covariance Matrices and Their Decompositions
After marginal transformation, NPDE computation estimates the covariance among simulated score vectors for a subject. A Cholesky factorization is efficient and convenient when the matrix is symmetric positive definite. If the estimated covariance is singular, Cholesky factorization fails. If it is nearly singular, rounding error can make the factorization fail or produce scores that are highly sensitive to small changes in the simulations.
Singularity has a structural meaning: at least one linear combination of the score components has zero simulated variance. This can occur when two measurements are exact duplicates, when the number of independent simulated vectors is too small relative to the number of observations, or when the model and design make an observation a deterministic combination of others. Near-singularity means some direction has very little simulated variation relative to other directions. Inverting that covariance can greatly magnify small discrepancies in that direction.
Worked Example: Detecting an Exactly Singular Covariance
Suppose three transformed score components have the estimated covariance \( C=\begin{pmatrix}1&1&0\\1&1&0\\0&0&1\end{pmatrix} \). The first two components are identical in the simulations: their variances are 1 and their covariance is also 1. The third is uncorrelated with them.
The first two rows of \(C\) are identical, so its determinant is 0 and it has no inverse. Equivalently, the vector \(v=(1,-1,0)^\mathsf{T}\) satisfies \(Cv=0\). Thus the simulated variance of the difference between the first two score components is zero. A Cholesky factorization requiring a positive definite matrix cannot be completed, and ordinary covariance whitening is undefined in that direction.
The first remedy is to check why the components are redundant: for example, whether the same observation was included twice or whether the data structure truly contains duplicate measurements. If they are duplicate records, remove the duplication and recompute. If the observations are genuinely distinct but deterministic under the simulated model, report that the covariance does not support a full-rank decorrelation; do not present a routine inverse as if it existed.
Adding a ridge term illustrates a possible numerical workaround, but also changes the calculation. With \(C_\lambda=C+0.01I\), the eigenvalues become 2.01, 1.01, and 0.01. The condition number, the largest eigenvalue divided by the smallest, is \(2.01/0.01=201\). The matrix is now invertible, but the variance assigned to the formerly zero-variance direction is an imposed value, not one learned from the simulations. The resulting NPDEs depend on that choice.
Near-Singularity and Sensitivity
A successful decomposition is not by itself proof of numerical reliability. A matrix can be positive definite but so poorly conditioned that small changes in estimated covariance lead to large changes in whitened scores. The condition number provides one warning: for a symmetric positive definite matrix it is the ratio of its largest to smallest eigenvalue. A large ratio signals that some directions are much less well determined than others. There is no universal cutoff independent of scale, precision, and purpose; inspect eigenvalues, warnings, and sensitivity rather than relying on one threshold.
Worked Example: A Small Eigenvalue Amplifies a Discrepancy
Consider the diagonal covariance matrix \(C=\operatorname{diag}(1,0.0011,0.4)\) and the marginal score vector \(e=(0.2,0.1414,-0.1)^\mathsf{T}\). Because the matrix is diagonal, its Cholesky factor has diagonal entries \(\sqrt{1}\), \(\sqrt{0.0011}\), and \(\sqrt{0.4}\). Whitening divides each component of \(e\) by the corresponding factor.
For the second component, \( \sqrt{0.0011}\approx 0.0331662 \), so \( 0.1414/\sqrt{0.0011}\approx 0.1414/0.0331662=4.263 \). The other components are \(0.2/\sqrt{1}=0.2\) and \(-0.1/\sqrt{0.4}\approx -0.1581\). Thus the decorrelated vector is approximately \((0.2,4.263,-0.1581)^\mathsf{T}\).
The covariance condition number is \(1/0.0011\approx 909.1\). The large second NPDE is driven by the small variance in that direction: a modest marginal discrepancy is magnified. It may reflect a genuinely unusual discrepancy, but it may also be sensitive to simulation error in the estimated covariance. Repeating the calculation with more simulations and inspecting the small-eigenvalue direction are therefore important before interpreting 4.263 as a stable diagnostic.
For illustration, adding \(0.01I\) changes the second variance from 0.0011 to 0.0111. The corresponding component becomes \(0.1414/\sqrt{0.0111}\approx 1.342\), rather than 4.263. This sensitivity shows why ridge regularization is not a neutral repair: it changes the implied scale of the discrepancy and can materially change the diagnostic.
Stable Computation and Practical Diagnostics
Use a decomposition suited to the matrix and retain its diagnostic information. For a well-conditioned positive definite covariance \(C=LL^\mathsf{T}\), solve the triangular system \(Lz=e\); do not explicitly calculate \(C^{-1}\) and multiply by \(e\). Solving the system is generally more stable and avoids an unnecessary matrix inverse. Verify symmetry up to numerical tolerance and check that the diagonal of the factor is positive and finite.
When Cholesky fails or the matrix appears ill-conditioned, an eigenvalue decomposition or singular-value decomposition can reveal the problematic directions. These methods help diagnose numerical rank, but they do not make a genuinely unidentified direction informative. If a small eigenvalue is truncated, the calculation effectively projects out that direction; if it is replaced with a floor, the resulting score depends on the chosen floor. State the rule and assess sensitivity.
Count boundary ranks and ties, confirm the plotting-position convention, and check that all normal scores are finite. Increase the number of simulations when rank resolution is inadequate for the intended diagnostic.
Check matrix dimensions, symmetry, diagonal variances, duplicate observations, eigenvalues, and condition number. Confirm that simulated vectors preserve the subject’s observation times and within-subject dependence.
Use Cholesky for a well-conditioned positive definite matrix. Use eigenvalue or singular-value diagnostics to investigate rank deficiency; do not silently treat a failed factorization as an ordinary NPDE result.
Compare results across adequate simulation sizes and any justified numerical tolerances. Record regularization or truncation choices, affected subjects, and whether conclusions change.
A further diagnostic is dependence on the order of observations. Cholesky whitening is tied to the chosen coordinate order: reordering components changes the triangular factor and may change individual whitened coordinates, even though the transformed covariance remains the identity in exact arithmetic. Use a fixed, documented order, typically chronological sampling order, when comparing results. A symmetric inverse-square-root transformation provides a different whitening convention, but it too must be specified; it does not resolve singularity or weakly estimated directions.
Remedies, Alternatives, and Limits
| Problem | Useful response | Important limitation |
|---|---|---|
| Boundary ranks or coarse percentiles | Increase the number of predictive simulations; use a documented finite-sample plotting position. | A finite score is still limited by the available simulation resolution. |
| Exact duplicate or redundant measurements | Correct duplicated records or define a nonredundant observation vector. | Do not remove genuinely distinct data merely to force a decomposition to succeed. |
| Ill-conditioned covariance | Increase simulation size, inspect small-eigenvalue directions, and compare stable decompositions. | More simulations may not fix structural redundancy or model-implied determinism. |
| Singular covariance with a justified computational need | Consider a documented eigenvalue truncation or ridge sensitivity analysis. | These choices alter the transformation; results are no longer the unmodified full-rank NPDE calculation. |
Increasing simulation size is often preferable when the difficulty arises from Monte Carlo error: it can improve percentile resolution and covariance estimation. It will not make an exactly redundant observation vector full rank. Dropping redundant components can be appropriate when their redundancy is real and understood. Regularization may be useful for a declared sensitivity analysis, but should not be used to conceal a failed assumption or to make an inconvenient score look acceptable.
Numerical warnings should be interpreted alongside the design and the predictive model. Check whether the issue recurs for particular subjects, sampling patterns, or outcome values, and whether it persists with more simulations. A numerical failure is not itself a goodness-of-fit result. Conversely, an extreme but stable NPDE should not be dismissed as numerical merely because it is inconvenient. The distinction rests on reproducibility, rank resolution, and the behavior of the estimated covariance.
Common Pitfall
A common mistake is to see an infinite score or failed Cholesky decomposition, add a small constant automatically, and then interpret the resulting NPDEs as if nothing changed. The reasoning is that the constant is “only numerical.” The correction is that a boundary-rank adjustment sets an empirical percentile convention, while a ridge changes covariance eigenvalues and therefore changes the scale and direction of whitening. Either can be defensible when specified and checked, but neither is an invisible repair. Record the choice, examine sensitivity, and distinguish the modified diagnostic from the unmodified calculation.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- With 199 simulations, an observation is below every simulated value. Using \( (r+t/2+0.5)/(M+1) \), with \(r=t=0\), what percentile and approximate normal score result?
- Why can a covariance matrix with a successful Cholesky factorization still yield unstable NPDEs?
- What does an exact zero eigenvalue mean for covariance whitening?
- Does adding \(0.01I\) leave the original NPDE calculation unchanged? Explain.
- Why should observation order be held fixed when comparing Cholesky-whitened NPDEs?
Solutions. (1) The percentile is \(0.5/(199+1)=0.0025\), giving \(\Phi^{-1}(0.0025)\approx -2.807\). (2) Cholesky success establishes positive definiteness at the computed precision, not good conditioning; a very small eigenvalue can amplify small errors. (3) At least one nonzero direction has zero variance, so the covariance has no ordinary inverse and full-rank whitening is undefined. (4) No. Adding \(0.01I\) changes the eigenvalues and the whitening scale, so the resulting scores depend on the ridge. (5) Cholesky whitening depends on coordinate order; changing order can change individual transformed components even when the joint covariance remains identity in exact arithmetic.