Cfrtjryk

This question is an extension of a related question about adding small probabilities. Suppose you have log-probabilities $ell_1 geqslant ell_2$, where the corresponding probabilities $exp(ell_1)$ and $exp(ell_2)$ are too small to be distinguished from zero in the initial computational facility being used (e.g., base R). We want to find the log-diffference of these probabilities, which we denote by:

$$ell_- equiv ln big( exp(ell_1) - exp(ell_2) big)$$

Questions: How can you effectively compute this log-difference? Can this be done in the base R? If not, what is the simplest way to do it with package extensions?

To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:

$$begin{equation} begin{aligned}
exp(ell_1) - exp(ell_2)
&= exp(ell_1) (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$

This result converts the difference to a product, which allows us to present the log-difference as:

$$begin{equation} begin{aligned}
ell_-
&= ln big( exp(ell_1) - exp(ell_2) big) \[6pt]
&= ln big( exp(ell_1) (1 - exp(-(ell_1 - ell_2))) big) \[6pt]
&= ell_1 + ln (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$

In the case where $ell_1 = ell_2$ we obtain the expression $ell_+ = ell_1 + ln 0 = -infty$. Using the Maclaurin series expansion for $ln(1-x)$ we obtain the formula:

$$begin{equation} begin{aligned}
ell_-
&= ell_1 - sum_{k=1}^infty frac{exp(-k(ell_1 - ell_2))}{k} quad quad quad text{for } ell_1 neq ell_2. \[6pt]
end{aligned} end{equation}$$

Since $exp(-(ell_1 - ell_2)) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $ell_1 - ell_2$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.

Implementation in base R: It is possible to compute this log-difference accurately in base R using the log1p function. This is a primitive function in the base package that computes the value of $ln(1+x)$ for an argument $x$ (with accurate computation even for $x ll 1$). This primitive function can be used to give a simple function for the log-difference:

logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }

Implementation with VGAM package: Machler (2012) analyses accuracy issues in evaluating the function $ln(1-exp(-|x|))$, and suggests that use of the base R functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp function in the VGAM package. This gives you the an alternative function for the log-difference:

logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }

The following workaround is often very useful for these sorts of problems:

1. subtract the smaller of l1 and l2 from each of l1 and l2 (effectively, we have multiplied the probabilities by some constant z = exp(min(l1,l2)))


2. Now compute the sum using standard functions (you can use log1p and expm1 if you want).


3. Afterwards, add the quantity you subtracted in step 1.

A simple example:

l1 <- -2000 ## exp(-2000) is computational zero

l2 <- -2002

z <- min(l1,l2)

l1 <- l1 - z

l2 <- l2 - z

## now we are guaranteed that one of them is zero

y <- log(exp(l1) - exp(l2))

y + z

returns the correct value -2000.145, which is equal to -2000 + log(1-exp(-2)).