在 R 中处理非常小的数字

Question

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

(0.1)^1000、0.2^(1200)，

然后将它们标准化，以便它们的总和为 1 即

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 这样的数据类型 可以让你得到比双精度更好的值。

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

Answer 1

从数学上来说，这些数字之一将是appx。零，还有另一个。你们的数字之间的差异很大，所以我什至想知道这是否有意义。

但一般来说，要做到这一点，您可以使用 R 底层的 logspace_add C 函数的想法。可以定义 logxpy ( =log(x+y) )logxpy ( =log(x+y) ) 当 lx = 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 减去 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 包处理非常大或小数字，本质上将乔里斯的答案包装成一种方便的形式。

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)

Answer 3

尝试任意精度包：

• Rmpfr "R MPFR -可靠的多精度浮点”
• Ryacas “R 接口“Yacas”计算机代数系统” - 也可能能够做到任意精度。

Answer 4

也许你可以将 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)))