r - 接近 1 的数的幂

Question

我猜有一些我无法找到的标准技巧：无论如何，我想以数值稳定的方式计算一个非常接近 1 的数字的大幂（想想 1-p，其中 p<1e-17） . 1-p 在我的系统上被截断为 1。

使用对数的泰勒展开，我得到以下界限

我有什么更聪明的办法吗？

Answer 1

您可以使用该功能log(1+x)更准确地计算。|x| <= 1log1p

一个例子：

> p <- 1e-17
> log(1-p)
[1] 0
> log1p(-p)
[1] -1e-17

还有一个：

> print((1+1e-17)^100, digits=22)
[1] 1
> print(exp(100*log1p(-1e-17)), digits=22)
[1] 0.9999999999999990007993

然而，在这里，我们受限于double基于类型的 FP 算术的准确性（请参阅What Every Computer Scientist Should Know About Floating-Point Arithmetic）。

另一种方法是使用例如Rmpfr（又名多精度浮点可靠）包：

> options(digits=22)
> library(Rmpfr)
> .N <- function(.) mpfr(., precBits = 200) # see the package's vignette
> (1-.N(1e-20))^100
1 'mpfr' number of precision  200   bits 
[1] 0.99999999999999999900000000000000005534172854579042829381053529

该包使用gslandmpfr库来实现任意精度的 FP 操作（当然，以较慢的计算速度为代价）。