Tutorials › Pharmacometric Workflows and Reproducible Modeling › Debugging Rounding Errors in a Likelihood Calculation

Debugging Model Failures · Tutorial 169 of 399

Debugging Rounding Errors in a Likelihood Calculation

Learn to locate precision loss in a likelihood calculation and verify whether a stable reformulation preserves the intended result.

Advanced 11 min read

What You'll Learn

  • Distinguish rounded input values, rounded intermediate calculations, and rounded display output.
  • Trace the effect of probability rounding on a Bernoulli log-likelihood contribution.
  • Reformulate sums of tiny likelihood terms with the log-sum-exp identity.
  • Recognize cancellation when computing small log-likelihood changes.
  • Check numerical results against precision, scale, and an independently computed formulation.

Why Rounding Can Change a Likelihood

A likelihood calculation can be mathematically correct and still give a misleading result if precision is lost along the way. This matters in pharmacometrics because likelihoods may combine many contributions, very small probabilities, and numerical predictions from a model. Rounding an output for a table is usually harmless; rounding an intermediate probability or log-likelihood contribution can change the calculation itself.

A useful starting point is to separate three operations. Input rounding replaces a value supplied to the calculation with a less precise value. Intermediate rounding rounds a value while calculating, such as a probability or a contribution to the objective. Display rounding changes only the printed representation after the calculation. Only the last of these is normally just a presentation choice. To debug a discrepancy, trace where each value was rounded and compare the value used in the calculation with the value printed.

Definition: For observations \(y_1,\ldots,y_n\) with likelihood contributions \(L_i(\theta)\), the likelihood is \(L(\theta)=\prod_{i=1}^{n}L_i(\theta)\), and the log-likelihood is \(\ell(\theta)=\sum_{i=1}^{n}\log L_i(\theta)\). A numerical rounding error occurs when finite-precision representation or an explicit rounding operation makes the computed value differ from the value intended by the mathematical expression. Whether that difference matters depends on its size relative to the scale and use of the result.

Log-likelihoods are generally preferable to multiplying many probabilities directly: a product of small positive numbers can underflow to zero, while its logarithm can remain representable. But using logs does not make every calculation stable. A probability may already have rounded to zero or one before its logarithm is taken, and subtracting nearly equal log-likelihoods can lose meaningful digits. The key is to preserve precision in the quantities that feed the log-likelihood and choose formulas that avoid unstable intermediate operations.

Worked Example: A Rounded Probability Becomes a Zero Likelihood

Worked Example: A Rounded Probability Becomes a Zero Likelihood

Consider one hypothetical Bernoulli observation \(y=0\). Let the model probability of \(y=1\) be \(p=0.9999999\), so the likelihood contribution for the observed zero is \(1-p=0.0000001\). The Bernoulli log-likelihood contribution is

$$ \ell(p\mid y=0)=\log(1-p). $$

Using the unrounded probability gives

$$ \ell=\log(1-0.9999999)=\log(0.0000001)\approx -16.1181. $$

Suppose instead that an intermediate step rounds \(p\) to six decimal places. Then \(p\) becomes \(1.000000\), so the calculation becomes \(\log(1-1)=\log(0)\), which is \(-\infty\). The rounded calculation has changed a finite log-likelihood contribution into a nonfinite one. This is not a small difference caused by formatting the final answer; the rounded probability has changed the value used in the likelihood.

A stable approach is to calculate the small probability directly when the model parameterization makes it available. For example, if \(q=1-p=10^{-7}\) is computed directly, the contribution is \(\log(q)=\log(10^{-7})\approx -16.1181\), without subtracting a rounded value of \(p\) from one. If \(p\) is retained with sufficient precision, \(\log1p(-p)\) is also a suitable way to calculate \(\log(1-p)\), particularly when \(p\) is close to zero. It cannot, however, recover information after \(p\) has already been rounded to exactly one.

The calculation is valid for a Bernoulli observation with a model probability strictly between zero and one. In this example, the assumed probability is in that range and the observed value is zero; the failure comes from rounding, not from an impossible observation under the stated model.

Stable Forms for Small Probabilities and Sums

Sometimes a model supplies a probability through a parameterization rather than as a separately calculated probability and complement. For a Bernoulli model with logit \(\eta=\log\{p/(1-p)\}\), the probability of a zero can be written as \(1-p=1/(1+\exp(\eta))\). Its log is

$$ \log(1-p)=-\log(1+\exp(\eta)). $$

For large positive \(\eta\), calculating \(p\) first and then subtracting it from one risks losing the small complement. A stable expression uses the softplus function, defined by \(\operatorname{softplus}(x)=\max(x,0)+\log1p(\exp(-|x|))\). Thus \(\log(1-p)=-\operatorname{softplus}(\eta)\). For \(\eta=\log(9999999)\), this gives \(-\log(10000000)\approx -16.1181\), retaining the finite log contribution without first rounding \(p\) to one.

Another common difficulty occurs when several positive likelihood terms are added and each term is tiny. If \(a_i\) are the log values of those terms, then the log of their sum is not the sum of their logs. It is

$$ \log\left(\sum_{i=1}^{m}\exp(a_i)\right) =a_{\max}+\log\left(\sum_{i=1}^{m}\exp(a_i-a_{\max})\right), \qquad a_{\max}=\max_i a_i. $$

This is the log-sum-exp identity. It preserves the result by shifting the exponentials so their largest exponent is zero rather than a large negative number. For finite \(a_i\), at least one shifted term is \(\exp(0)=1\), so the sum inside the logarithm is at least one and does not underflow to zero.

Justification: Since \(a_{\max}\) is common to every term, \(\sum_i\exp(a_i)=\sum_i\exp(a_{\max})\exp(a_i-a_{\max})=\exp(a_{\max})\sum_i\exp(a_i-a_{\max})\). Taking the logarithm and using \(\log(uv)=\log(u)+\log(v)\) gives the identity. The shift changes the arithmetic form, not the mathematical value.

Worked Example: Adding Likelihood Terms on the Log Scale

Worked Example: Adding Likelihood Terms on the Log Scale

Suppose a hypothetical calculation requires the log of the sum of two positive terms, whose log values are \(a_1=-1000\) and \(a_2=-1001\). Directly exponentiating these values in ordinary floating-point arithmetic can underflow: both \(\exp(-1000)\) and \(\exp(-1001)\) may be represented as zero. Taking the log of that computed sum then yields \(\log(0)=-\infty\).

Apply log-sum-exp with \(a_{\max}=-1000\):

$$ \begin{aligned} \log\{\exp(-1000)+\exp(-1001)\} &=-1000+\log\{\exp(0)+\exp(-1)\}\\ &=-1000+\log(1+0.3678794)\\ &\approx -1000+0.3132617\\ &\approx -999.6867. \end{aligned} $$

The stable result is finite even though the unscaled terms are too small for direct exponentiation in the assumed arithmetic. This example assumes that the two log terms themselves are available as finite numbers; log-sum-exp does not repair an earlier error that made either input log value incorrect or nonfinite. It also does not decide whether the terms represent an appropriate likelihood model: it only evaluates their sum more reliably.

When Nearly Equal Values Are Subtracted

Rounding can also matter when the target is a small change rather than a large log-likelihood value. If two fitted log-likelihoods are close, subtracting them may discard leading digits that agree. Likewise, a formula such as \(\log(1+\delta)\) can lose a small \(\delta\) when \(1+\delta\) rounds to one. A numerically stable function such as \(\log1p(\delta)\) evaluates the intended small change directly.

Worked Example: Preserving a Small Log-Likelihood Change

Suppose a hypothetical likelihood ratio is \(1+\delta\), with \(\delta=10^{-17}\). The associated log-likelihood difference is \(\log(1+\delta)\). In arithmetic where \(1+10^{-17}\) rounds to \(1\), the direct calculation returns \(\log(1)=0\). The small positive change has disappeared.

Using \(\log1p\) instead gives

$$ \log1p(10^{-17})\approx 10^{-17}. $$

The approximation follows from the expansion \(\log(1+\delta)=\delta-\delta^2/2+\delta^3/3-\cdots\) for \(|\delta|<1\). With \(\delta=10^{-17}\), the second term has magnitude \(5\times10^{-35}\), so the leading term is \(10^{-17}\) to the displayed precision. The direct formula loses the difference because it forms \(1+\delta\) first; the stable function is designed to preserve the small increment.

This example illustrates a numerical mechanism, not a claim that every likelihood ratio this close to one is scientifically meaningful. Whether a tiny computed change matters depends on the model, the data, and the decision being made. The numerical calculation should be checked before interpreting the change.

Tracing and Verifying a Rounding Problem

When a likelihood result changes after rounding, inspect the calculation in sequence rather than changing model settings at random. Record the values before and after each operation, including the precision of inputs and intermediate values. Then compare the original formula with a mathematically equivalent stable form. Agreement between independently evaluated forms is evidence that the discrepancy has been localized; it is not proof that the likelihood model itself is correct.

1
Identify the first differing value.
Compare the unrounded calculation with the rounded one at each intermediate step. Check especially probabilities near zero or one, very small exponentials, and differences between similar values.
2
Separate calculation from display.
Determine whether the low-precision number was merely printed or was actually passed into a later calculation. Reformatting a stored high-precision result should not change the underlying likelihood.
3
Use a stable equivalent expression.
Use direct log probabilities or \(\log1p\) for complements and small increments, and log-sum-exp for sums of exponentials. Preserve the unrounded inputs to those expressions.
4
Check sensitivity and reproducibility.
Repeat the calculation with additional retained digits or a higher-precision implementation if available. Confirm that the conclusion is not an artifact of the chosen precision, and report enough digits to show the meaningful difference without implying unsupported accuracy.

Diagnostics include an abrupt change from a finite value to \(-\infty\), a probability printed as zero or one when a finite log contribution is expected, or a result that changes substantially when only an intermediate rounding setting changes. A discrepancy that remains when stable formulas and full-precision inputs are used may instead arise earlier in the model calculation, such as from an inaccurate numerical prediction or an invalid likelihood contribution. Those are different failure sources and should not be mislabelled as rounding.

Choosing an Approach and Knowing Its Limits

ApproachUseful whenMain limitation
Direct probability arithmeticProbabilities and sums are comfortably within the representable range and are not near cancellation.Products can underflow; complements can round to zero; small changes can disappear.
Log probabilities and log-sum-expLikelihoods involve products, very small probabilities, or sums of positive terms.Inputs must still be computed accurately; the transformation does not correct model or integration errors.
Higher precisionA stable reformulation is insufficient or an independent precision check is needed.It may cost more computation and does not guarantee a correct formula or scientifically appropriate model.
Rounded outputPresenting a result at a readable precision after calculation is complete.It is not a remedy for rounded intermediate values and may conceal meaningful differences if too few digits are shown.

A stable reformulation is usually preferable to simply demanding more digits: it avoids a known source of numerical loss while retaining ordinary precision for the rest of the calculation. Higher precision can serve as a useful check, but agreement at higher precision cannot establish that the likelihood has been specified correctly. Conversely, a tiny difference in the last printed digit is not automatically a failure. Assess its size relative to the intended comparison and the precision justified by the inputs and computation.

Common Mistake

Common mistake: “The probability printed as 1.000000, so the observed zero has likelihood zero; the model must be impossible.” The displayed number may be a rounded representation of a value below one. Inspect the unrounded value and the calculation path. If rounding occurred before subtraction or logging, calculate the small complement directly or use a stable log-probability expression. A displayed value alone cannot establish that the underlying likelihood contribution is zero.

The reverse mistake is also possible: assuming that a stable formula makes every result trustworthy. Stable arithmetic can evaluate the stated expression accurately, but it cannot repair incorrect data, an inappropriate probability model, or an inaccurate upstream prediction. Numerical verification and model checking answer related but distinct questions.

Check Your Understanding

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

  1. A Bernoulli observation is \(y=0\) and its model probability is \(p=0.9999\). Write its log-likelihood contribution and give an approximate value.
  2. Why can \(\log1p(-p)\) fail to help if \(p\) has already been rounded to exactly one?
  3. Use log-sum-exp to evaluate \(\log\{\exp(-20)+\exp(-21)\}\).
  4. What is the difference between display rounding and intermediate rounding?

Solutions. For question 1, \(\ell=\log(1-0.9999)=\log(0.0001)\approx -9.2103\). For question 2, once the stored value is exactly one, its complement is zero; a stable function cannot reconstruct the discarded information. The complement must be preserved or computed another way before that loss. For question 3, the result is \(-20+\log(1+\exp(-1))\approx -20+0.3132617=-19.6867\). For question 4, display rounding changes only how a completed value is printed; intermediate rounding changes a value used by subsequent operations and can therefore alter the calculated likelihood.

Clinical Biostats