我猜测有一些标准的技巧我没能找到:无论如何,我想以数值稳定的方式计算一个非常接近1(比如1-p,其中p < 1e-17)的大数次幂。在我的系统上,1-p被截断为1。
使用对数的泰勒展开式,我得到以下界限:
exp(-np-np^2/2) ≤ exp(n*log(1-p)) ≤ exp(-np)
是否有更聪明的方法?
使用对数的泰勒展开式,我得到以下界限:
exp(-np-np^2/2) ≤ exp(n*log(1-p)) ≤ exp(-np)
是否有更聪明的方法?
通过使用 log1p
函数,您可以更准确地计算 |x| <= 1
的 log(1+x)
。
例如:
> 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算术准确性的限制(请参见计算机科学家应知道的浮点算术)。
另一种方法是使用例如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
gsl
和mpfr
库来实现任意精度的浮点数运算(当然会导致计算速度变慢)。