Resta de probabilidades muy pequeñas - ¿Cómo se calcula?

Question

Esta pregunta es una extensión de una pregunta relacionada con la adición de pequeñas probabilidades . Supongamos que se tienen probabilidades logarítmicas $\ell_1 \geqslant \ell_2$ donde las probabilidades correspondientes $\exp(\ell_1)$ y $\exp(\ell_2)$ son demasiado pequeños para distinguirse de cero en la instalación computacional inicial que se utiliza (por ejemplo, la base R ). Queremos encontrar la diferencia logarítmica de estas probabilidades, que denotamos por:

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

Preguntas: ¿Cómo se puede calcular efectivamente esta diferencia logarítmica? ¿Puede hacerse en la base R ? Si no es así, ¿cuál es la forma más sencilla de hacerlo con extensiones de paquetes?

Answer 1

Para ver cómo tratar este tipo de diferencias, primero observamos un resultado matemático útil sobre las diferencias de exponenciales:

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

Este resultado convierte la diferencia en un producto, lo que nos permite presentar la diferencia logarítmica como

$$\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}$$

En el caso de que $\ell_1 = \ell_2$ obtenemos la expresión $\ell_+ = \ell_1 + \ln 0 = -\infty$ . Utilizando el Ampliación de la serie Maclaurin para $\ln(1-x)$ obtenemos la fórmula:

$$\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}$$

Desde $\exp(-(\ell_1 - \ell_2)) < 1$ los términos de esta expansión disminuyen rápidamente (más rápido que el decaimiento exponencial). Si $\ell_1 - \ell_2$ es grande, entonces los términos disminuyen particularmente rápido. En cualquier caso, esta expresión nos permite calcular la suma logarítmica a cualquier nivel de precisión deseado, truncando la suma infinita a un número deseado de términos.

---

Implementación en R base: Es posible calcular con precisión esta diferencia logarítmica en base R utilizando el log1p función. Es una función primitiva del paquete base que calcula el valor de $\ln(1+x)$ para un argumento $x$ (con un cálculo preciso incluso para $x \ll 1$ ). Esta función primitiva puede utilizarse para obtener una función simple para la diferencia logarítmica:

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

Aplicación con VGAM paquete: Machler (2012) analiza los problemas de precisión en la evaluación de la función $\ln(1-\exp(-|x|))$ y sugiere que el uso de la base R funciones puede implicar una pérdida de precisión. Es posible calcular esta diferencia logarítmica con mayor precisión al utilizar el log1mexp en la función VGAM paquete. Esto le da una función alternativa para la diferencia logarítmica:

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

Answer 2

La siguiente solución suele ser muy útil para este tipo de problemas:

1. restar el menor de l1 y l2 de cada uno de l1 y l2 (efectivamente, hemos multiplicado las probabilidades por alguna constante z = exp(min(l1,l2)) )
2. Ahora calcule la suma utilizando funciones estándar (puede utilizar log1p y expm1 si quieres).
3. Después, añade la cantidad que has restado en el paso 1.

Un ejemplo sencillo:

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

devuelve el valor correcto -2000.145 que es igual a -2000 + log(1-exp(-2)) .