compute probabilities from log scores

Question

I'm using R and have the following task. Let's say I have access to log(a) and log(b), and I want to compute a/(a+b). Further suppose that log(a) and log(b) are very negative, e.g. -1000, then how do I compute a/(a+b) in the most numerically efficient way? Keep in mind that due to underflow, exp(log(a)) will output a 0.

How should I do this in R?

Consider the following test case:

log_a = log(2e-100)
log_b = log(4e-100) 

I want a function which can take these as log scores and output 0.33.

Answer 1

The following functions are transcripts of the code in this SO post by user Joris Meys, hence the name of function joris.

naive <- function(lx, ly) {
  ex <- exp(lx)
  ey <- exp(ly)
  ex/(ex + ey)
}

logxpy <- function(lx, ly) {
  max(lx, ly) + log1p(exp(-abs(lx - ly)))
}
joris <- function(lx, ly) {
  exp(lx - logxpy(lx, ly))
}

log_a <- log(2e-100)
log_b <- log(4e-100)
log_x <- -1000
log_y <- -1001

naive(log_a, log_b)
#> [1] 0.3333333
joris(log_a, log_b)
#> [1] 0.3333333

naive(log_x, log_y)
#> [1] NaN
joris(log_x, log_y)
#> [1] 0.7310586

^{Created on 2022-05-03 by the reprex package (v2.0.1)}