我猜想有一些我找不到的标准技巧:无论如何,我想以数值稳定的方式计算一个非常接近1的大幂(认为1-p,其中p
使用对数的泰勒展开式,我得到以下边界
有什么我可以做的更聪明的吗?
最佳答案
您可以使用log(1+x)
函数为|x| <= 1
更准确地计算log1p
。
一个例子:
> 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
该程序包使用
gsl
和mpfr
库来实现任意精度的FP操作(当然会以较慢的计算速度为代价)。