Why Equivalent Formulas Can Behave Differently
In pharmacometric calculations, a probability or likelihood contribution may be expressed as the difference between two values that are each close to one another. The exact difference can be meaningful, but floating-point arithmetic stores only a finite number of significant digits. If rounding has already obscured the small difference between the operands, subtracting them cannot recover it. This is cancellation error.
The problem is not that subtraction is inherently inaccurate. Subtracting two exactly represented, well-separated numbers can be reliable. The risk is that the operands are close in size, often because an earlier calculation rounded each one. Difference expressions appear in pharmacometric work, for example, when calculating the probability that an event time lies within a narrow interval or evaluating a model term involving two similar exponentials. An algebraic rearrangement can preserve the small quantity directly rather than obtain it by subtracting two rounded, nearly equal values.
A useful first check is the relative size of a difference. For nonzero \(x\) and \(y\), define its subtraction condition number as
When \(x\) and \(y\) are close, the denominator is small and \(\kappa_{\mathrm{sub}}\) is large. Small relative errors in the operands can then become much larger relative errors in \(x-y\). This indicator describes sensitivity to operand error; it does not, by itself, prove that a particular software calculation is inaccurate.
A Floating-Point View of Cancellation
In the standard floating-point model, a basic operation on ordinary, representable values is often described as \(\operatorname{fl}(a\mathbin{\circ} b)=(a\mathbin{\circ} b)(1+\delta)\), where \(\circ\) is an arithmetic operation and the relative rounding error \(\delta\) has magnitude at most the unit roundoff \(u\). This description assumes the result is in the normal range and that no exceptional behavior occurs. It is a model for understanding rounding, not a guarantee for every library function.
Suppose the intended operands are \(x\) and \(y\), but their computed values have small errors. The resulting error in their difference is influenced by errors on the scale of \(|x|+|y|\), whereas the desired answer has scale \(|x-y|\). When \(x\) and \(y\) nearly agree, their ratio can be large. For example, errors around the last few digits of each operand may be negligible relative to either operand but substantial relative to their difference.
A stable reformulation avoids creating that vulnerable subtraction, when possible. The functions expm1 and log1p are designed for this purpose: expm1(\(x\)) evaluates \(e^x-1\), and log1p(\(x\)) evaluates \(\log(1+x)\), with better accuracy when \(x\) is near zero than a direct evaluation of the displayed subtraction or sum followed by a standard function. Stable does not mean exact; it means the calculation is arranged to retain accuracy in the situation at hand.
Worked Example: Evaluating a Small Exponential Difference
Consider the dimensionless quantity \(q=1-e^{-z}\) for \(z=10^{-16}\). This form can occur when computing a small event probability over a short interval. Compare the direct calculation with an expm1 calculation, assuming binary64 arithmetic and a correctly rounded exponential function.
The exact value is \(1-e^{-10^{-16}}\), approximately \(9.9999999999999995\times10^{-17}\). In binary64, \(e^{-10^{-16}}\) rounds to the representable number immediately below 1, namely \(1-2^{-53}\). The direct subtraction therefore returns
Its relative error is approximately \((1.1102230-1.0000000)/1.0000000=0.110223\), or 11.0%. The result is not zero. For comparison, at \(z=10^{-17}\), the correctly rounded binary64 value of \(e^{-z}\) is 1, so direct subtraction returns zero even though the exact \(1-e^{-z}\) is positive.
Instead, use \(q=-\operatorname{expm1}(-z)\). The expm1 function evaluates the small difference without first rounding \(e^{-z}\) to a value nearly equal to 1. At \(z=10^{-16}\), the result is approximately \(10^{-16}\), retaining the meaningful size of the quantity. This example illustrates the numerical principle; actual results depend on the quality and precision of the software library.
Stable Evaluation of Likelihood Contributions
Suppose a hypothetical event time \(T\) follows an exponential distribution with rate \(\lambda>0\), measured in \(\mathrm{h}^{-1}\). Its survival function is \(P(T>t)=e^{-\lambda t}\) for \(t\geq0\). For \(0\leq a<b\), the probability that the event occurs in the interval \((a,b]\) is the difference
When \(b-a\) is small, both terms can be close and direct subtraction may lose relative accuracy. Factoring out the first exponential gives an equivalent expression in real arithmetic:
The factorization isolates the small interval width in the expm1 term. If the likelihood requires the logarithm of this probability, the same algebra gives a direct log contribution: \(-\lambda a+\log[-\operatorname{expm1}(-\lambda(b-a))]\). Thus, one need not first calculate the difference of two similar survival probabilities and then take its logarithm.
Worked Example: An Exponential Event-Time Interval
For a clearly hypothetical calculation, let \(\lambda=0.200\ \mathrm{h}^{-1}\), \(a=50.000000\ \mathrm{h}\), and \(b=50.000006\ \mathrm{h}\). The model and its assumptions are those just stated: a constant positive exponential rate and an interval with \(0\leq a<b\). The interval width is \(6.0\times10^{-6}\ \mathrm{h}\), so \(z=\lambda(b-a)=1.2\times10^{-6}\).
The direct probability formula subtracts \(e^{-10}\) and \(e^{-10.0000012}\), two nearby values. In the rearranged formula, evaluate the small factor using expm1:
This is a dimensionless probability. The corresponding log probability, useful when this interval probability is a likelihood contribution, is
The interval ordering and positive rate ensure this probability is positive, so its logarithm is defined. The stable form improves the calculation of the small difference; it does not change the assumed exponential model or correct misspecified interval data.
Difference of Exponentials and a Limiting Case
A second recurring pattern is a quotient of a difference of exponentials by a difference in rates. For distinct positive rates \(k_a\) and \(k_e\), both in \(\mathrm{h}^{-1}\), consider the time-dependent expression
Let \(\Delta=k_a-k_e\) and \(z=\Delta t\). Factoring \(e^{-k_e t}\) and then using expm1 yields
The final ratio has a removable singularity at \(z=0\). Define \(\phi(z)=-\operatorname{expm1}(-z)/z\) for \(z\ne0\), and define \(\phi(0)=1\), its continuous limit. The stable form \(e^{-k_e t}t\phi(z)\) then also gives the limiting value when \(z=0\), including \(t=0\). This qualification matters: the displayed quotient with division by \(z\) is not itself defined at \(z=0\). The continuous extension is a separate definition.
Worked Example: Close Rate Constants
Consider the hypothetical rates \(k_e=0.1000000\ \mathrm{h}^{-1}\), \(k_a=0.1000001\ \mathrm{h}^{-1}\), and time \(t=12.0\ \mathrm{h}\). These rates are positive and distinct, as required for the original quotient. Here \(\Delta=1.0\times10^{-7}\ \mathrm{h}^{-1}\) and \(z=\Delta t=1.2\times10^{-6}\).
The stable calculation is
The units are hours: the exponential and ratio are dimensionless, and the explicit factor \(t\) supplies the unit. Directly evaluating the original numerator subtracts two exponentials near \(e^{-1.2}\), then divides by a small rate difference. The rewritten form computes the small relative change through expm1 instead. It is algebraically equivalent for the stated distinct rates, but it is generally better conditioned in floating-point arithmetic for this close-rate case.
Using Logarithms of Differences
Some likelihood calculations need the logarithm of a positive difference \(e^A-e^B\), where \(A>B\). Directly forming that difference and then taking its logarithm can discard precision when \(A\) and \(B\) are close. Factoring out \(e^A\) gives
Since \(B-A<0\), the inner quantity is positive. For \(d=B-A\) near zero, evaluate \(\log(1-e^d)\) as \(\log[-\operatorname{expm1}(d)]\). When \(d\) is sufficiently negative, a library function for \(\log(1+x)\), such as log1p, can evaluate the expression as \(\operatorname{log1p}(-e^d)\). These are algebraic alternatives selected for different numerical regimes; the required condition \(A>B\) must still hold. If the ordering fails, the real logarithm of the difference is not defined as written.
Diagnostics, Alternatives, and Limitations
A suspicious difference should prompt a numerical check, not an immediate change to the pharmacometric model. Warning signs include a small probability becoming zero or negative, a result changing materially under an algebraically equivalent formulation, or a likelihood contribution changing when the calculation is repeated at higher precision. Check the signs and ordering of the operands, the units and scale of the inputs, and whether a value is truly close to zero or merely rounded for display. Compare the direct expression with a stable reformulation at representative parameter values, including the limiting cases relevant to the model.
The principal alternatives have different costs and uses:
| Approach | Useful when | Main caution |
|---|---|---|
| Direct subtraction | Operands are not close, or the difference is not sensitive to small operand errors. | Can lose relative accuracy for nearly equal operands. |
| Algebraic reformulation using expm1 or log1p | The expression has a known structure such as \(1-e^x\) or \(\log(1+x)\), especially near zero. | Requires correct algebra, domain conditions, and treatment of limiting cases. |
| Higher precision arithmetic | A stable reformulation is unavailable or is being used as a numerical comparison. | Can increase computational cost and does not fix incorrect inputs or a flawed model. |
These methods address arithmetic precision, not every source of error. A stable expression cannot recover information lost when input values were rounded before evaluation. It also cannot resolve whether the model or observed data are appropriate. Special-function implementations vary in quality, and a formula that is stable for a small argument may need a different evaluation strategy for other argument ranges. The practical requirement is to validate the chosen computation over the parameter and data ranges used by the analysis.
Identify whether the operands are close relative to their magnitudes, and check whether the resulting difference is a small probability or likelihood factor.
Verify signs, ordering, parameter restrictions, and whether a denominator can be zero. State a continuous limiting value only when it exists and is justified.
Use expm1, log1p, factorization, or higher precision as appropriate; compare results and test important limiting cases.
Check Your Understanding
Work through each question, then compare with the solutions below.
- Why can direct subtraction of two nearly equal floating-point quantities have poor relative accuracy, even if each operand is individually accurate?
- For a positive, small \(z\), give a stable evaluation of \(1-e^{-z}\). What is the value of the direct binary64 subtraction at \(z=10^{-17}\), assuming correctly rounded exponential evaluation?
- For an exponential event-time model with rate \(\lambda>0\), write a stable expression for \(P(a<T\leq b)\) when \(0\leq a<b\).
- In the rate-difference expression, why must the form involving \(-\operatorname{expm1}(-z)/z\) specify what happens at \(z=0\)?
- Does agreement between a direct and a reformulated calculation establish that the statistical model is correct? Explain.
Solutions. 1. Operand errors are on the scale of the operands, but the answer is their much smaller difference; those errors can therefore be large relative to the answer. 2. Use \(-\operatorname{expm1}(-z)\). At \(z=10^{-17}\), the rounded exponential is 1, so direct subtraction returns 0 even though the exact result is positive. 3. Use \(e^{-\lambda a}[-\operatorname{expm1}(-\lambda(b-a))]\); the stated rate and interval conditions make the probability positive. 4. The ratio divides by \(z\), so it is undefined at zero. Its continuous value is 1, but that extension must be defined explicitly. 5. No. Numerical agreement supports the reliability of the evaluated expression for those inputs; it does not establish that the model assumptions or data are correct.