r - 1に近い数の累乗

Question

私が見つけることができなかったいくつかの標準的なトリックがあると思います: とにかく、数値的に安定した方法で 1 に非常に近い数 (p<1e-17 である 1-p を考えてください) の大きな累乗を計算したいのです。 . 私のシステムでは、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 演算の精度に制限があります (すべてのコンピューター科学者が浮動小数点演算について知っておくべきことを参照してください)。

もう 1 つの方法は、たとえばRmpfr(別名 Multiple Precision Floating-Point Reliable) パッケージを使用することです。

> 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 演算を実装します (もちろん、計算速度は遅くなります)。