r语言 - 从日志分数计算概率

Question

我使用R，并有以下任务。假设我有log(a)和log(b)，我想计算a/(a+b)。进一步假设log(a)和log(b)是非常负的，例如-1000，那么我如何以最有效的数字方式计算a/(a+b)?请记住，由于下溢，exp(log(a))将输出一个0。

在R中我该怎么做呢?

考虑以下测试用例:

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

我想要一个函数，可以把这些作为日志分数和输出0.33。

Answer 1

以下函数是用户Joris Meys在这篇SO帖子中代码的转录本，因此函数名称为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

^{在2022-05-03由reprex包(v2.0.1)创建}