Subtracting very small probabilities - How to compute? 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?
 A: The following workaround is often very useful for these sorts of problems:


*

*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)))

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

*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)).
A: 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); }

