23

RIで、私が説明できない奇妙な行動を見つけています。ここの誰かが説明できることを望んでいます。100の価値があると思います！この大きな数です。

予想される動作を示すコンソールからの数行...

>factorial( 10 )
[1] 3628800
>prod( 1:10 )
[1] 3628800
> prod( as.double(1:10) )
[1] 3628800
> cumprod( 1:10 )
[1]       1       2       6      24     120     720    5040   40320  362880 3628800

しかし、100を試してみると！私は次のようになります（結果の数値が約14桁でどのように異なり始めるかに注意してください）：

> options(scipen=200) #set so the whole number shows in the output
> factorial(100)
[1] 93326215443942248650123855988187884417589065162466533279019703073787172439798159584162769794613566466294295348586598751018383869128892469242002299597101203456
> prod(1:100)
[1] 93326215443944102188325606108575267240944254854960571509166910400407995064242937148632694030450512898042989296944474898258737204311236641477561877016501813248
> prod( as.double(1:100) )
[1] 93326215443944150965646704795953882578400970373184098831012889540582227238570431295066113089288327277825849664006524270554535976289719382852181865895959724032
> all.equal( prod(1:100) , factorial(100) , prod( as.double(1:100) ) )
[1] TRUE

「既知の」数100に設定された変数に対していくつかのテストを実行すると！次に、次のように表示されます。

# This is (as far as I know) the 'true' value of 100!
> n<- as.double(93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000)
> factorial(100) - n
[1] -1902315522848807765998160811905210717565551993186466795054798772271710903343294674760811531554315419925519536152107160826913610179566298858520576
> prod(1:100) - n
[1] -48777321098687378615337456715518223527321845979140174232174327494146433419058837814379782860367062049372295798771978482741374619988879457910784
> prod(as.double(1:100)) - n
[1] 0

最終結果はゼロと評価されますが、値100！に設定された変数がどこにあるprod( as.double( 1:100 ) )かを正しく評価しても、返される数値は期待どおりに表示されません。prod( as.double( 1:100 ) ) - nn

誰かが私にこの振る舞いを説明できますか？私がx64システムを使用しているので、私が知る限り、オーバーフローなどに関連するべきではありません。以下のバージョンとマシン情報：

> .Machine$double.xmax
[1] 1.798e+308
> str( R.Version() )
List of 14
 $ platform      : chr "x86_64-apple-darwin9.8.0"
 $ arch          : chr "x86_64"
 $ os            : chr "darwin9.8.0"
 $ system        : chr "x86_64, darwin9.8.0"
 $ status        : chr ""
 $ major         : chr "2"
 $ minor         : chr "15.2"
 $ year          : chr "2012"
 $ month         : chr "10"
 $ day           : chr "26"
 $ svn rev       : chr "61015"
 $ language      : chr "R"
 $ version.string: chr "R version 2.15.2 (2012-10-26)"
 $ nickname      : chr "Trick or Treat"

誰かが私にこれを説明できますか？Rがすべてを正しく実行することは間違いありません。これは、おそらくuseRに関連しています。prod( as.double( 1:100 ) ) - n私が気になっていることを正しく評価しているので、プロジェクトオイラー 問題20を実行しているので、正しい数字を表示する必要があることを指摘するかもしれません。

ありがとう

4

4 に答える 4

15

これは、の最大値ではdoubleなく、その精度と関係があります。

100!158の有効な（10進数の）桁があります。IEEE double（64ビット）には仮数用に52ビットのストレージスペースがあるため、小数点以下約16桁の精度を超えると丸め誤差が発生します。

ちなみに、100!あなたが疑ったように、実際には、

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

したがって、計算されたRの値はすべて正しくありません。

Rはわかりませんが、比較する前all.equal()にこれら3つの値すべてをsに変換しているようです。そのため、それらの違いは失われます。float

于 2013-01-14T11:00:40.910 に答える
13

でのテストでは、期待どおりの結果が得られall.equalません。2つの値all.equalしか比較できません。3番目の引数は、位置的にに一致します。これにより、比較操作の許容範囲が与えられます。あなたへのあなたの呼び出しにおいて、それの許容範囲を与えてください、それは間違いなく、ばかげて異なる値に対して真である比較につながります：toleranceall.equal100!

> all.equal( 0, 1000000000, prod(as.double(1:100)) )
[1] TRUE

しかし、あなたがそれに2つの引数だけを与えたとしても、例えば

all.equal( prod(1:100), factorial(100) )

TRUEデフォルトの許容値はであるため、それでも生成されます.Machine$double.eps ^ 0.5。たとえば、2つのオペランドは約8桁に一致する必要がありますが、これは間違いなく当てはまります。一方、許容値をに設定すると0、3つの可能な組み合わせのどちらも比較から等しくなりません。

> all.equal( prod(1:100), factorial(100), tolerance=0.0 )
[1] "Mean relative difference: 1.986085e-14"
> all.equal( prod(1:100), prod( as.double(1:100) ), tolerance=0.0 )
[1] "Mean relative difference: 5.22654e-16"
> all.equal( prod(as.double(1:100)), factorial(100), tolerance=0.0 )
[1] "Mean relative difference: 2.038351e-14"

また、Rに200の有効数字を印刷するように指示したからといって、それらがすべて正しいことを意味するわけではないことにも注意してください。実際、1 / 2 ^ 53の10進数は約53桁ですが、意味があると見なされるのは最初の16桁だけです。

これにより、「真の」値との比較にも欠陥が生じます。これを観察してください。Rが提供する最後の桁は次のfactorial(100)とおりです。

...01203456

あなたはそれから引きますn、ここnで100の「本当の」値です！したがって、最後に24個のゼロが必要であり、したがって、差も同じ桁で終了する必要がfactorial(100)あります。しかし、それは次のように終わります。

...58520576

これは、これらすべての数字が重要ではないことを示しているだけであり、実際にそれらの値を調べるべきではありません。

100を正確に表すには、525ビットの2進精度が必要です。-これはの10倍の精度ですdouble。

于 2013-01-14T13:24:29.603 に答える
7

発生している動作をグラフィカルに説明するために、3番目の回答を追加します。基本的に、階乗計算の倍精度は22！までで十分です。その後、実際の値からますます発散し始めます。

50！付近では、factorial（x）とprod（1：x）の2つのメソッドがさらに区別され、後者は、ご指摘のとおり、「実際の」ファクターにより近い値を生成します。

Rの階乗計算精度

添付コード：

# Precision of factorial calculation (very important for the Fisher's Exact Test)
library(gmp)
perfectprecision<-list()
singleprecision<-c()
doubleprecision<-c()
for (x in 1:100){
    perfectprecision[x][[1]]<-factorialZ(x)
    singleprecision<-c(singleprecision,factorial(x))
    doubleprecision<-c(doubleprecision,prod(1:x))
}


plot(0,col="white",xlim=c(1,100),ylim=c(0,log10(abs(doubleprecision[100]-singleprecision[100])+1)),
        ,ylab="Log10 Absolute Difference from Big Integer",xlab="x!")
for(x in 1:100) {
    points(x,log10(abs(perfectprecision[x][[1]]-singleprecision[x])+1),pch=16,col="blue")
    points(x,log10(abs(perfectprecision[x][[1]]-doubleprecision[x])+1),pch=20,col="red")
}
legend("topleft",col=c("blue","red"),legend=c("factorial(x)","prod(1:x)"),pch=c(16,20))
于 2013-07-11T23:29:46.783 に答える
1

さて、あなたはそれが呼び出すことを本体から知ることができますfactorial、それはgammaを呼び出します.Primitive("gamma")。どのように.Primitive("gamma")見えますか？このように。

大きな入力の場合、.Primitive("gamma")の動作はそのコードの198行目にあります。それは呼んでいます

exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
            ((2*y == (int)2*y)? stirlerr(y) : lgammacor(y)));

これは単なる概算です。

ちなみに、その例としての使用に関する記事。Rmpfrfactorialしたがって、問題を解決しようとしている場合は、「Rmpfrライブラリを使用するだけ」です。

于 2015-09-03T17:06:04.317 に答える