如何处理两个以上超小概率之和的对数

问题描述 投票：0回答：2

出于某种邪恶的原因，我需要计算 500 个超小概率之和的对数，每一项由

计算 
dmvnorm(X[,i], mean=rep(0,3), sigma=diag(3))

有时上面的代码会由于下溢而返回 0，但使用对数就可以了。

dmvnorm(X[,i], mean=rep(0,3), sigma=diag(3), log=TRUE)

我知道我可以用数学方法处理两个术语：

log(p1 + p2) = log(p2) + log(1 + p1/p2)
。但是我们可以将这种方法推广到更多术语吗？有谁对此有更多经验吗？

顺便说一句，我写了一个递归函数来计算这个。从数学上来说，它是有效的。但我认为这不实用。

MESSY <- function (pv) 
{
  N <- length(pv)
  if (N==1)
    {return(pv[1])}
  else
    {w <- pv[N]
     pv <- pv[1:N-1]-w
     return(w + log(1 + exp(MESSY(pv))))
    }
}

因为有些 p 非常小，我只能使用 

w=log(p)
，所以我们有 
log(exp(w1)+exp(w2)) = w2 + log(1+exp(w1-w2))
和 
log(exp(w1)+exp(w2)+exp(w3)) = w3 + log(1 + exp(w1-w3) + exp(w2-w3))
两个和三个术语。

r probability logarithm underflow
2个回答
2
投票

这个函数是从R源代码中的内部

logspace_add
函数翻译过来的这里

logspace_add <- function(logx,logy) { 
    pmax(logx,logy) + log1p(exp(-abs(logx - logy))) 
}

不一定是最有效的，但您应该能够使用 

Reduce()
:

对 >2 个元素执行此操作 
logspace_add_mult <- function(x) {
    Reduce(logspace_add, x)
}

快速测试（使用足够大不会下溢的值，以便我们可以比较常规计算和对数空间计算的结果）。

x <- c(1e-4,1e-5,1e-6)
logspace_add_mult(log(x))
## [1] -9.10598
log(sum(x))
## [1] -9.10598

据我所知，这或多或少是日志空间添加的标准方法。使用除此本土实现之外的其他实现的优点是 (1) 代码成熟度和测试以及 (2) 速度（至少对于 

logspace_add_mult
版本；我怀疑原生 C（或不管怎样）实施
logspace_add
）。 Brobdingnag 包 使用类似的表示：

library(Brobdingnag)
brob(log(x))
## [1] +exp(-9.2103) +exp(-11.513) +exp(-13.816)
sum(brob(log(x)))
## [1] +exp(-9.106)
log(as.numeric(sum(brob(log(x)))))
## [1] -9.10598

在 Python 中，numpy 有 logaddexp，但这只能成对使用：您可以使用 

functools.reduce()
来概括它，如上所述。

import numpy as np
import functools as f
x = np.array([1e-4,1e-5,1e-6])
f.reduce(np.logaddexp,np.log(x))

这可能比 R 中的 

Reduce()
快一点。

0
投票

添加到@BenBolker的答案 A method to vectorize 

logspace_add
:

logspace_sum1 <- function(logx) {
  while (length(logx) > 1L) {
    if (length(logx)%%2) {
      m <- matrix(logx[-1], ncol = 2, byrow = TRUE)
      logx <- c(logx[1], pmax(m[,1], m[,2]) + log1p(exp(-abs(m[,1] - m[,2]))))
    } else {
      m <- matrix(logx, ncol = 2, byrow = TRUE)
      logx <- pmax(m[,1], m[,2]) + log1p(exp(-abs(m[,1] - m[,2])))
    }
  }
  
  logx
}

另一种选择是减去 

log(sum(exp()))
之前向量中的最大值，然后将其加回来：

logspace_sum2 <- function(logx) {
  if (max(logx) < log(.Machine$double.xmin)) {
    m <- max(logx)
    log(sum(exp(logx - m))) + m
  } else log(sum(exp(logx)))
}

测试并计时我们的选择：

library(Rmpfr)
library(Brobdingnag)

funs <- c("log(sum(exp(logx)))",
          "Reduce(logspace_add, logx)",
          "logspace_sum1(logx)",
          "logspace_sum2(logx)",
          "as.numeric(log(sum(exp(as.brob(logx)))))",
          "as.numeric(log(sum(exp(mpfr(logx, 128)))))"
)

第一个向量将因天真的

log(sum(exp()))
而下溢。

set.seed(1688732183)

logx <- Rfast::dmvnorm(matrix(runif(3e4, 35, 40), 1e4, 3), numeric(3), diag(3), TRUE)
summary(logx)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   -2386   -2179   -2114   -2115   -2049   -1854

创建一个矩阵来保存结果。

logsumx <- matrix(0, length(funs), 1, dimnames = list(funs, NULL))

microbenchmark::microbenchmark(
  list = setNames(lapply(lapply(1:length(funs), \(f) parse(text = sprintf("logsumx[%d] <- %s", f, funs[f]))), \(x) bquote({..(x)}, splice = TRUE)), funs),
  times = 1
)
#> Unit: microseconds
#>                                        expr      min       lq     mean   median       uq      max neval
#>                         log(sum(exp(logx)))     44.5     44.5     44.5     44.5     44.5     44.5     1
#>                  Reduce(logspace_add, logx) 609559.8 609559.8 609559.8 609559.8 609559.8 609559.8     1
#>                         logspace_sum1(logx)  11191.8  11191.8  11191.8  11191.8  11191.8  11191.8     1
#>                         logspace_sum2(logx)    343.2    343.2    343.2    343.2    343.2    343.2     1
#>    as.numeric(log(sum(exp(as.brob(logx)))))   7448.5   7448.5   7448.5   7448.5   7448.5   7448.5     1
#>  as.numeric(log(sum(exp(mpfr(logx, 128)))))  83300.0  83300.0  83300.0  83300.0  83300.0  83300.0     1

print(logsumx, digits = 15)
#>                                                         [,1]
#> log(sum(exp(logx)))                                     -Inf
#> Reduce(logspace_add, logx)                 -1854.05150072123
#> logspace_sum1(logx)                        -1854.05150072123
#> logspace_sum2(logx)                        -1854.05150072123
#> as.numeric(log(sum(exp(as.brob(logx)))))   -1854.05150072123
#> as.numeric(log(sum(exp(mpfr(logx, 128))))) -1854.05150072123

第二个数据集不会下溢：

logx <- Rfast::dmvnorm(matrix(runif(3e4, 20, 25), 1e4, 3), numeric(3), diag(3), TRUE)
summary(logx)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  -932.4  -804.2  -763.4  -764.8  -725.0  -608.9

microbenchmark::microbenchmark(
  list = setNames(lapply(lapply(1:length(funs), \(f) parse(text = sprintf("logsumx[%d] <- %s", f, funs[f]))), \(x) bquote({..(x)}, splice = TRUE)), funs),
  times = 1
)
#> Unit: microseconds
#>                                        expr      min       lq     mean   median       uq      max neval
#>                         log(sum(exp(logx)))    582.7    582.7    582.7    582.7    582.7    582.7     1
#>                  Reduce(logspace_add, logx) 641831.3 641831.3 641831.3 641831.3 641831.3 641831.3     1
#>                         logspace_sum1(logx)   1759.7   1759.7   1759.7   1759.7   1759.7   1759.7     1
#>                         logspace_sum2(logx)   4088.9   4088.9   4088.9   4088.9   4088.9   4088.9     1
#>    as.numeric(log(sum(exp(as.brob(logx)))))   2671.1   2671.1   2671.1   2671.1   2671.1   2671.1     1
#>  as.numeric(log(sum(exp(mpfr(logx, 128)))))  72447.6  72447.6  72447.6  72447.6  72447.6  72447.6     1

logsumx[] <- sapply(funs, \(x) eval(parse(text = x)))
print(logsumx, digits = 15)
#>                                                         [,1]
#> log(sum(exp(logx)))                        -608.329353240961
#> Reduce(logspace_add, logx)                 -608.329353240961
#> logspace_sum1(logx)                        -608.329353240961
#> logspace_sum2(logx)                        -608.329353240961
#> as.numeric(log(sum(exp(as.brob(logx)))))   -608.329353240961
#> as.numeric(log(sum(exp(mpfr(logx, 128))))) -608.329353240961

在不溢出的方法中，

logspace_sum2
的性能最好。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.