18

我需要计算一个非常小的数字列表,例如

(0.1)^1000, 0.2^(1200),

然后将它们归一化,以便它们总和为一个,即

a1 = 0.1^1000, a2 = 0.2^1200

我想计算 a1' = a1/(a1+a2), a2'=a2(a1+a2)。

当我得到 a1=0 时,我遇到了下溢问题。我怎样才能解决这个问题?理论上我可以处理日志,然后 log(a1) = 1000*log(0.l) 将是一种表示 a1 而没有下溢问题的方法 - 但为了规范化我需要得到 log(a1+a2) -我无法计算,因为我不能直接表示 a1 。

我正在用 R 编程——据我所知,c# 中没有像 Decimal 这样的数据类型,它可以让你获得比双精度值更好的值。

任何建议将不胜感激,谢谢

4

4 回答 4

16

从数学上讲,其中一个数字将是 appx。零,另一个。您的数字之间的差异很大,所以我什至想知道这是否有意义。

但是一般来说,要做到这一点,您可以使用logspace_addR 引擎盖下的 C 函数中的想法。可以定义logxpy ( =log(x+y) )whenlx = log(x)ly = log(y)as :

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

这意味着我们可以使用:

> la1 <- 1000*log(0.1)
> la2 <- 1200*log(0.2)

> exp(la1 - logxpy(la1,la2))
[1] 5.807714e-162

> exp(la2 - logxpy(la1,la2))
[1] 1

如果您有更多数字,也可以递归调用此函数。请注意,1 仍然是 1,而不是 1 minus 5.807...e-162。如果您确实需要更高的精度并且您的平台支持 long double 类型,您可以使用例如 C 或 C++ 编写所有代码,然后稍后返回结果。但如果我是对的,R 可以——目前——只处理正常的双打,所以当结果显示时,你最终会再次失去精度。


编辑 :

为你做数学:

log(x+y) = log(exp(lx)+exp(ly))
         = log( exp(lx) * (1 + exp(ly-lx) )
         = lx + log ( 1 + exp(ly - lx)  )

现在你只取最大的 lx,然后你得到表达式logxpy()

编辑2:为什么要取最大值呢?很简单,确保您在 exp(lx-ly) 中使用负数。如果 lx-ly 变得太大,则 exp(lx-ly) 将返回 Inf。这不是一个正确的结果。exp(ly-lx) 将返回 0,从而获得更好的结果:

说 lx=1 和 ly=1000,然后:

> 1+log1p(exp(1000-1))
[1] Inf
> 1000+log1p(exp(1-1000))
[1] 1000
于 2011-04-27T11:40:03.857 回答
8

Brobdingnag软件包处理非常大或非常小的数字,基本上将 Joris 的答案包装成一种方便的形式。

a1 <- as.brob(0.1)^1000
a2 <- as.brob(0.2)^1200
a1_dash <- a1 / (a1 + a2)
a2_dash <- a2 / (a1 + a2)
as.numeric(a1_dash)
as.numeric(a2_dash)
于 2011-04-27T12:59:14.443 回答
2

尝试任意精度包:

  • Rmpfr“R MPFR - 多精度浮点可靠”
  • Ryacas“R Interface to the 'Yacas' Computer Algebra System” - 也可以做任意精度。
于 2011-04-27T12:29:40.103 回答
1

也许您可以将 a1 和 a2 视为分数。在您的示例中,使用

a1 = (a1num/a1denom)^1000  # 1/10
a2 = (a2num/a2denom)^1200  # 1/5

你会到达

a1' = (a1num^1000 * a2denom^1200)/(a1num^1000 * a2denom^1200 + a1denom^1000 * a2num^1200)
a2' = (a1denom^1000 * a2num^1200)/(a1num^1000 * a2denom^1200 + a1denom^1000 * a2num^1200)

可以使用 gmp 包计算:

library(gmp)
a1 <- as.double(pow.bigz(5,1200) / (pow.bigz(5,1200)+ pow.bigz(10,1000)))
于 2011-04-27T11:49:57.053 回答