处理R中的非常小的数字

Question

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

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

然后将它们归一化,使它们总结为一个ie

a1 = 0.1 ^ 1000,a2 = 0.2 ^ 1200

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

我遇到了下溢问题,因为我得到a1 = 0.我怎么能绕过这个？从理论上讲,我可以处理日志,然后log(a1)= 1000*log(0.l)将是一种表示没有下溢问题的方法 - 但为了规范化我需要得到log(a1 + a2) - 我无法计算,因为我不能直接代表a1.

我用R编程 - 据我所知,c#中没有数据类型如Decimal,它可以让你比双精度值更好.

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

Answer 1

在数学上说,其中一个数字将是appx.零,另一个.你的数字之间的差异是巨大的,所以我甚至想知道这是否有意义.

但是为了做到这一点,你可以使用来自logspace_addR的引擎盖下的C函数的想法.可以定义logxpy ( =log(x+y) )何时lx = log(x)和ly = log(y)作为:

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减去5.807...e-162.如果您真的需要更高的精度并且您的平台支持长双精度类型,您可以使用例如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

Answer 2

该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)