Stable Computation of Censored Likelihoods

Why Stable Computation Matters

A below-quantification-limit (BQL) record does not report an exact concentration. It reports that the concentration fell below a known limit. Under a specified residual-error model, its likelihood contribution is therefore a probability. For an extreme lower-tail event, that probability can be so small that a computer represents it as zero, even though its logarithm is finite. A zero contribution is not a harmless approximation: its logarithm is negative infinity and can derail likelihood evaluation and optimization.

The central computational rule is to evaluate the log cumulative probability directly, rather than first calculating a tiny probability and then taking its logarithm. This tutorial focuses on that calculation at a fixed set of model parameters and predictions. A correct censored likelihood still depends on choosing the appropriate error distribution, limit, and conditioning structure.

Definition: If a record is known only to satisfy \(Y<L\), where \(L\) is its lower limit of quantification, its censored likelihood contribution under model parameters \(\theta\) is \(P(Y<L\mid\theta)\). Its log-likelihood contribution is \(\log P(Y<L\mid\theta)\). For a continuous distribution, using \(Y\leq L\) gives the same probability.

The calculation assumes that the residual distribution used in the model is appropriate, that the relevant limit \(L\) is known and correctly assigned, and that the BQL status corresponds to the event \(Y<L\). If a likelihood is formed as a product over records, the records must be independent conditional on the quantities being conditioned on, or their joint likelihood must otherwise be modelled correctly. Stable arithmetic cannot repair a misspecified distribution, an incorrect limit, or an inappropriate independence assumption.

From a BQL Event to a Log CDF

For additive normal residual error, suppose that at the evaluation point the concentration \(Y\) has distribution \(Y\sim N(F,\sigma^2)\), where \(F\) is the model-predicted mean concentration and \(\sigma>0\) is the residual standard deviation. Standardizing the event gives

P(Y<L)=P\left(\frac{Y-F}{\sigma}<\frac{L-F}{\sigma}\right) =\Phi(z),\qquad z=\frac{L-F}{\sigma},

where \(\Phi\) is the standard normal cumulative distribution function. Thus the required contribution is \(\log\Phi(z)\). When \(F\) is far above \(L\), \(z\) is strongly negative and \(\Phi(z)\) is very small.

For a lognormal residual model, suppose \(\log Y\sim N(\log F,\sigma^2)\), with \(F>0\) the median concentration and \(\sigma\) the standard deviation on the log scale. For \(L>0\),

P(Y<L)=\Phi\left(\frac{\log L-\log F}{\sigma}\right).

This too is a normal log-CDF calculation, but its standardized argument is obtained on the log-concentration scale. Using the additive-normal argument in a lognormal model, or vice versa, evaluates the wrong likelihood.

Worked Example: An Extreme Additive-Normal Tail

Suppose a hypothetical BQL record has \(L=0.010\) mg/L, with \(F=0.410\) mg/L and additive residual standard deviation \(\sigma=0.0100\) mg/L. The normal model gives

z=\frac{0.010-0.410}{0.0100}=-40.0,\qquad \ell_{\mathrm{BQL}}=\log\Phi(-40.0).

Directly evaluating \(\Phi(-40)\) in ordinary floating-point arithmetic commonly returns zero. Instead use the lower-tail expansion described below. With \(x=40\), its leading terms give

\log\Phi(-x)\approx -\frac{x^2}{2}-\log x-\frac{1}{2}\log(2\pi) +\log\left(1-\frac{1}{x^2}+\frac{3}{x^4}-\frac{15}{x^6}+\cdots\right).

At \(x=40\), the bracketed correction begins \(1-1/1600+3/2{,}560{,}000-15/4{,}096{,}000{,}000\), approximately \(0.999376\). Therefore,

\log\Phi(-40)\approx -800-\log(40)-0.918939+\log(0.999376) \approx -804.608.

The log contribution is finite, even though exponentiating it produces a probability too small for ordinary double-precision representation. This calculation is appropriate only if the additive-normal model and the stated values are the intended ones.

Stable Ways to Evaluate the Normal Log CDF

A robust implementation should call a numerical routine designed to return \(\log\Phi(z)\) directly (often named a normal log-CDF or log-normal-CDF routine). Such routines use different formulas in different regions to avoid underflow and loss of precision. When a suitable routine is unavailable, a useful exact identity for \(z<0\) uses the scaled complementary error function,

\operatorname{erfcx}(x)=e^{x^2}\operatorname{erfc}(x),\qquad \log\Phi(z)=-\frac{z^2}{2}+\log\operatorname{erfcx}\left(\frac{-z}{\sqrt{2}}\right)-\log 2.

This follows from \(\Phi(z)=\tfrac12\operatorname{erfc}(-z/\sqrt{2})\). The scaling in \(\operatorname{erfcx}\) compensates for the rapidly decaying exponential that makes direct evaluation of \(\operatorname{erfc}\) problematic in the tail. The formula is only useful if the implementation of \(\operatorname{erfcx}\) itself is numerically reliable over the required range.

An asymptotic expansion can also provide a check or a fallback. It is important to understand why it has the form shown in the first example. Write \(x>0\), and let \(\phi(t)=(2\pi)^{-1/2}e^{-t^2/2}\). Since \(\phi'(t)=-t\phi(t)\), integration by parts gives

\Phi(-x)=\int_x^\infty\phi(t)\,dt =\frac{\phi(x)}{x}-\int_x^\infty\frac{\phi(t)}{t^2}\,dt.

Applying the same integration-by-parts step to the remaining integral produces successively smaller inverse powers of \(x\), yielding

\Phi(-x)\sim\frac{\phi(x)}{x} \left(1-\frac{1}{x^2}+\frac{3}{x^4}-\frac{15}{x^6}+\cdots\right), \qquad x\to\infty.

Taking logarithms gives the log-tail expansion. This is an asymptotic, not generally convergent, series: adding terms indefinitely is not a reliable strategy. Truncate while the terms are decreasing, and prefer a tested log-CDF routine for production calculations. For moderately negative arguments, the exact scaled-error-function formula or a validated library routine is preferable to relying on a short tail expansion.

Worked Example: A Lognormal BQL Contribution

Consider a hypothetical record with \(L=0.0500\) mg/L, lognormal median prediction \(F=0.200\) mg/L, and log-scale residual standard deviation \(\sigma=0.300\). The model implies

z=\frac{\log(0.0500)-\log(0.200)}{0.300} =\frac{\log(0.25)}{0.300} \approx -4.621.

The contribution is \(\log\Phi(-4.621)\), not the log density at \(L\). A reliable log-CDF calculation gives approximately \(-13.169\), corresponding to a BQL probability of about \(1.91\times10^{-6}\). As a check, the tail expansion with \(x=4.621\) gives

-\frac{x^2}{2}-\log x-\frac12\log(2\pi) +\log\left(1-\frac{1}{x^2}+\frac{3}{x^4}-\frac{15}{x^6}+\cdots\right) \approx -13.169.

The expansion is less extreme here than in the previous example, so the agreement should be checked against a trusted log-CDF routine rather than treated as exact from just a few terms. The arithmetic uses the lognormal model's log-scale standard deviation and a positive quantification limit, as required.

Accumulate Likelihoods on the Log Scale

For conditionally independent BQL records with standardized arguments \(z_1,\ldots,z_m\), the joint likelihood contribution is a product of probabilities. Its log is a sum:

L_{\mathrm{BQL}}=\prod_{i=1}^{m}\Phi(z_i), \qquad \log L_{\mathrm{BQL}}=\sum_{i=1}^{m}\log\Phi(z_i).

The product can underflow even when each individual probability is representable. Summing log contributions avoids that failure and is also the standard way to combine likelihood contributions numerically. The conditional-independence assumption matters: if it does not hold, simply summing marginal log probabilities is not the correct joint log likelihood.

Worked Example: Many Small Contributions

Suppose 100 hypothetical, conditionally independent BQL records each have \(z=-10\) under a fixed model evaluation. A stable normal log-CDF gives \(\log\Phi(-10)\approx-53.2313\). The joint log contribution is

\log L_{\mathrm{BQL}} =100(-53.2313)=-5323.13.

The probability for one record is about \(7.620\times10^{-24}\). Multiplying 100 such probabilities produces a number on the order of \(10^{-2312}\), below the representable range of standard double-precision arithmetic; direct multiplication can therefore return zero. The sum of the 100 log probabilities remains a finite value, approximately \(-5323.13\). This avoids numerical underflow, provided each individual log-CDF is itself computed stably.

Gradients, Diagnostics, and Failure Modes

Likelihood maximization often requires derivatives as well as likelihood values. For the standard normal log-CDF,

\frac{d}{dz}\log\Phi(z)=\frac{\phi(z)}{\Phi(z)}.

This ratio is sometimes called the inverse Mills ratio. It is positive and grows in the far lower tail. If \(z\) depends on a model parameter \(\theta\), the chain rule gives

\frac{d}{d\theta}\log\Phi(z(\theta)) =\frac{\phi(z)}{\Phi(z)}\frac{dz}{d\theta}.

Do not calculate the ratio by separately evaluating \(\phi(z)\) and \(\Phi(z)\) when \(\Phi(z)\) may underflow. If a stable log-CDF is available, calculate it as \(\exp(\log\phi(z)-\log\Phi(z))\), or use a specialized stable derivative routine. For example, at \(z=-40\) the ratio is approximately \(40.025\), not infinity. A zero CDF from naive evaluation would instead create an undefined ratio or an invalid gradient.

Useful diagnostics include checking that every log contribution is finite for finite inputs, comparing implementations at representative tail values, and checking that the total log likelihood equals the sum of the individual log contributions. Also inspect gradients and optimizer messages. Abrupt jumps to negative infinity, zero probabilities at plausible inputs, NaN gradients, or estimates that change markedly when a numerical routine is replaced are signs to investigate—not evidence that the model has found a better fit.

Approach	What it does	Main concern
Directly calculate \(\log(\Phi(z))\)	Computes the CDF and then takes its logarithm.	Fails if the CDF underflows to zero; may also lose precision near one.
Use a dedicated log-CDF or scaled-error-function formula	Computes the log probability without first forming a tiny probability.	Requires a tested implementation with adequate tail accuracy.
Use an asymptotic tail expansion	Approximates the extreme lower tail from inverse powers of \(-z\).	Truncation matters; it is an approximation and should not be extended blindly.
Replace zero probabilities with a fixed floor	Prevents a numerical zero by imposing a minimum contribution.	Changes the likelihood and can introduce artificial flat regions or parameter-dependent bias.

For arguments close to zero or positive and large, a stable log-CDF routine is still preferable. In the upper tail, \(\Phi(z)\) can round to one, so taking its logarithm may return zero even when the log probability is slightly negative. A routine may handle this using a survival probability and a stable log-one-minus-probability calculation. This is a different precision problem from lower-tail underflow, but it is another reason not to assemble the result from a naively rounded CDF.

A practical implementation should preserve the model's actual censoring probability. Substituting a fixed small probability when a calculation returns zero may conceal a software problem and alter both the likelihood and its derivatives. Replacing BQL values with numbers is also not a stable way to compute the censored likelihood: it changes the observation being modelled. Stable evaluation addresses arithmetic; it does not resolve questions about the suitability of the residual model, the observed censoring process, or the information available in heavily censored data.

Key takeaway: For a BQL record, calculate \(\log P(Y<L)\) under the specified residual distribution directly. For normal and lognormal residual models this is a normal log-CDF; evaluate it with a stable routine, sum log contributions, and treat zeros or non-finite gradients as numerical failures to diagnose rather than values to patch.

Check Your Understanding

Answer each question, then compare with the solutions below.

Under additive normal residual error, what standardized value determines the BQL probability for prediction \(F\), standard deviation \(\sigma\), and limit \(L\)?
Why is \(\log\Phi(-40)\) preferable to calculating \(\Phi(-40)\) and then taking its logarithm?
For a lognormal model with \(F=0.100\) mg/L, \(L=0.0500\) mg/L, and \(\sigma=0.200\) on the log scale, calculate the standardized value \(z\). Is the event in the lower tail?
What is the joint log contribution for 20 conditionally independent BQL records, each with log contribution \(-12.0\)?
Why is replacing every computed zero probability with a fixed floor not an equivalent numerical method?

Solutions. (1) \(z=(L-F)/\sigma\), and the probability is \(\Phi(z)\). (2) The CDF at \(-40\) can underflow to zero, whereas its log probability is finite, approximately \(-804.608\). (3) \(z=\log(0.0500/0.100)/0.200=\log(0.5)/0.200\approx-3.466\), so it is a lower-tail event. (4) Add the log contributions: \(20(-12.0)=-240\). (5) A fixed floor changes the likelihood and its gradient; it can conceal underflow rather than evaluate the intended probability accurately.