| English | 日本語 |
この文書は 英語版 の訳です。内容が食い違うときは英語版が正です。
公開しているすべての関数を R 4.6.2 と照合してテストしています (R にない erf と erfc は mpmath と照合)。期待値は oracle/ の R スクリプトで生成し、test/fixtures/ に保存しています (1 ケースにつき JSON 1 行)。そのため、テストの実行に R は要りません。テストの入力には必ず、典型的な値、遠い裾 (p = 1e-300)、境界 (0, 1, df → ∞)、NaN と Infinity、負の値、巨大なパラメータを含めています。
substitute)、理由を件数とともに出力します。rStats.js が意図して R と異なるケースは、テストファイルで理由を書いたうえで除外し (knownDifference)、differences-from-r.md に載せています。テストの出力は、フィクスチャごとに、R と一致したケース、参照値と一致したケース、後退誤差で合格したケース、除外したケースの数を報告します。
許容誤差は上限であって、目標ではありません。R が正しい範囲では、rStats.js は R と同じビットを目指します。別の計算に切り替えるのは、ふだんの計算の結果が有限でない、範囲外、または桁落ちしているときだけです。そのため、遠い範囲の修正で通常の範囲の結果が変わることはありません。R が誤っている範囲では、その範囲に限って rStats.js は正しい値を返し、違いを differences-from-r.md に記録します。
比較の詳細は次のとおりです。
| 0 に近い対数密度と対数確率は | 差 | / max( | 期待値 | , 1) で測ります。対数スケールでは 0 付近の値に特別な意味がないためです。 |
test/fixtures/<name>.ref.json に保存します (npm run refresh-references)。各 oracle は write_fixture(..., ref = function(row) ...) で参照値の関数を登録します。pnbinom の 4 件は、Rmpfr の参照値が誤っているか得られないため、代わりに mpmath (60 桁、betainc) で計算しました: 両側の pnbinom(1e8, 2.5, mu = 1e6) (Rmpfr は p = size/(size + mu) を 256 ビットで作ってから 1 − (反対側の裾) を取るため、8e-106 という結果には雑音しか残りません)、pnbinom(300, 50, 0.3, lower.tail = FALSE)、pnbinom(30, 3, 1 - 1e-10, log.p = TRUE)。nbinom について npm run refresh-references をやり直すと、最初のものは誤った Rmpfr の値に戻ってしまいます。両方のパラメータが大きい不完全ベータ関数 (a と b が 15〜1e10 の pbeta と pbinom。一様漸近展開の領域とその周辺、両側の裾と対数) も mpmath と照合します (scripts/mpmath/largebeta.py、audit_pbeta_large と audit_pbinom_large、許容誤差 1e-14。最大誤差は 9.7e-16 と 3.1e-15)。そこでは mpmath の betainc が数分かかるので、参照値は連分数 DLMF 8.17.22 を 50 桁と 80 桁で計算したもので (両者が 1e-28 まで一致したものだけを使う)、数値積分とも照合しています。一方のパラメータが大きく他方が小さい場合 (a が 15〜1e10、b が 0.1〜40 の pbeta とその鏡像、df が 30〜1e12 の pt。大きな a に対する展開の領域とその周辺) は、小さいほうのパラメータの枠での正の項の級数を 60 桁と 90 桁で計算したものを参照値とします (scripts/mpmath/smallb.py、audit_pbeta_smallb と audit_pt_smallb、許容誤差 1e-14。最大誤差は 5.4e-15 と 1.2e-15 で、不完全ベータ関数の級数によるもの。展開を通るケースは 1.2e-15 以内)。その領域の R のフィクスチャ (oracle/largea.R) には、すべてのケースに Rmpfr の参照値があります。R が約 1.0e-14 ずれていて、大きな a に対する展開で rStats.js がより正確になった (4e-16 以内) 5 件には、R が手元になく Rmpfr の参照値を計算できなかったため、mpmath の参照値 (100 桁と 200 桁、betainc、引数の double の値から。scripts/mpmath/mismatch_ref.py) を付けました: pbeta(0.1, 2.5, 1000) (上側の裾と、下側の対数)、pnbinom(30, 1000, mu = 2, log.p = TRUE)、pf(3, 10, 100, ncp = 0) (上側の裾と、下側の対数)。非心分布では、ポアソン混合級数の和または数値積分を Rmpfr で計算します (oracle/_noncentral.R)。実用的な時間で参照値を計算できない領域 (df = 1000 または 1e5 の非心 t、df2 = 1e9 かつ x > 30 の非心 F) はグリッドから外しています。reference フィールドに保存します。高精度が必要なのは平均、分散、相関係数だけなので、これは速く済みます。t 検定では、分散と標準誤差を Rmpfr (256 ビット) で計算しますが、平均は先に double へ丸めます。R も rStats.js も、estimate に出力する丸めた平均の差を使うためです。そのうえで p 値と分位点は R の pt/qt から求めます。Fisher の r × c 検定はすべての表を列挙し、gmp で厳密に合計します。Kendall の正確検定は Rmpfr (1024 ビット) で反転数を数えます。choose、lchoose、factorial、lfactorial、digamma、trigamma、psigamma、ベッセル関数、dmultinom、および R にない gammainc、erfinv、erfcinv、perm。これら 4 つの R のフィクスチャは pgamma × gamma、pracma の定義、Rmpfr で計算します): R のフィクスチャに加えて mpmath とも照合します (scripts/mpmath/*.py、60〜800 桁、test/fixtures/<name>.mpmath.json)。許容誤差は同じ 1e-14 です (erfinv、erfcinv、perm は 1e-15)。mpmath の値は数学的な関数なので、この照合では rCompatible: false を使います。R が不正確な点や失敗する点の参照値 (<name>.ref.json) は、Rmpfr ではなく mpmath (100〜800 桁、引数は倍精度の値) で計算しました。項が打ち消し合う箇所 (x < 0 の digamma と psigamma、振動する領域の J と Y、負の非整数次数の I) では、打ち消し合う項の大きさ (フィクスチャの scale) に対する相対誤差で測ります。x < 0 の psigamma は 2e-14 まで許します。反射項が cot(pi x) の deriv + 1 次の多項式で、cot の誤差を最大 101 倍にするためです (R との照合でも同じ)。besselJ と besselY の R との照合では、答えが打ち消し合う項よりはるかに小さいケースは、mpmath の値との差が scale に対して 1e-14 以内なら合格とします (knownDifference)。そこでは相対誤差で見ると R のほうが近いためです。log1pmx、lgamma1p、log1mexp、log1pexp。R では C API にしかありません): フィクスチャは R の出力ではなく、Rmpfr で計算した正しく丸めた値 (512 ビット。x が非常に小さい lgamma1p では 1 + x が厳密になるようにビット数を増やす) で、許容誤差 1e-14 で比較します (観測された最大の誤差は lgamma1p の 1.7e-15)。lmer): lme4 の推定値は最適化が止まった点で、θ では最適点から相対でおおむね 1e-9〜1e-5 離れています。そのため fixture には、厳密な最適点での参照値を持たせています。oracle (oracle/lmer.R、lmermore.R、lmerboundary.R) は lme4 の θ から出発し、基準の勾配が 0 になる点を 128 ビットの演算 (Rmpfr。基準は論文の式を密行列で、勾配は中心差分で) のニュートン法で求めます。最適点が境界上にある要素は、そこでの片側の微分が正であることを確かめてから 0 に固定します。そのうえで、その θ ですべての結果の項目を計算します。与えた θ での逸脱度 (lmerDeviance) には 200 ビットの参照値があります。構造 (lmerTerms: 計画行列、Z、共分散の因子の配置) は lme4 と完全に一致させます。solve と行列式には、gmp の有理数で計算した厳密な参照値があります (どの double も有理数で正確に表せます)。分解はそうではなく、後退誤差と直交性で確かめます (下記参照)。| 分野 | 相対許容誤差 |
|---|---|
| 特殊関数、密度関数、分布関数 | 1e-14 (非心 t: 1e-13、ptukey: 2e-14) |
| base R にない分布 (逆ガンマ、Kumaraswamy、Pareto、三角、Laplace、逆正弦) | d/p/q は 1e-14。extraDistr (逆正弦は R の beta(1/2, 1/2) を (a, b) に移したもの) と mpmath に対して |
| 分位点関数 | 1e-14、または後退誤差が R 以下。離散分布の分位点は完全一致か定義による検証 |
| 乱数 | 離散値、sample、r2dtable、L’Ecuyer-CMRG の runif と各ケース後の状態 (.Random.seed) は完全一致。連続値は 1e-12 以内。消費した一様乱数の個数を厳密に確かめる |
| 記述統計 | 1e-14 (打ち消し合いから決まる 0 に近い歪度と尖度は絶対誤差で測る) |
| 仮説検定 | すべての数値フィールドで 1e-12 (cor: 1e-15、pAdjust: 完全一致。ansariTest の正規近似の区間など、許容誤差 1e-4 の uniroot で求める根: 1e-4)。テキストフィールドは完全一致 |
線形モデル (lm、confint、predict) |
すべての数値フィールドで 1e-12 |
混合モデル (lmer) |
厳密な最適点の参照値に対して、logLik、AIC、BIC、deviance、REMLcrit は 1e-12。θ、固定効果とその共分散と標準誤差、σ、VarCorr、ranef、当てはめ値と残差は、それぞれの最大の絶対値の 1e-9 (観測した最大は 1.2e-12。境界のすぐ内側の平らな最適点)。参照値で 0 の θ の要素はちょうど 0。lme4 に対しては、基準の値が lme4 より大きくないことだけ。与えた θ での逸脱度は 1e-12 |
| 線形代数 | 1e-14 (crossprod、norm、backsolve、chol、qr、svd による 2-ノルム、行列演算、seq、lu) または 1e-13 (solve、determinant、rcond、svd、kappa、ginv、eigen)。結果のノルムに対する相対値で測る。prcomp は 1e-12。ピボット、ランク、行列式の符号は完全一致 |
| 数値計算 | 1e-15 (approx)、1e-14 (integrate、optimize、uniroot、barylag、hermiteInterp、ode、および絶対誤差の下限 1e-12 付きの grad、jacobian、hessian)、1e-13 (spline、splinefun、splinefunH) |
許容誤差は分野ごとに実測をもとに決めました。R または参照値に対して観測された最大の誤差は、許容誤差をわずかに下回っています。
ptukey (2e-14): R の ptukey の精度は約 8 桁なので、ほぼすべてのケースを参照値と比べます。参照値は 160 ビットの台形公式で、刻み幅を半分にしたときの値の変化が 1e-18 未満だった場合だけ使います。誤差のほとんどは 6e-15 未満で、最大は ptukey(0.5, 1000, 10) ≈ 1.3e-129 での 1.7e-14 です。ここでは log W の丸め (ε |
log W | ) が支配的です。 |
Math.log、Math.exp、Math.pow、Math.tan の最後の数ビットが V8 と glibc の libm で異なるため、連続値は完全には一致しません。観測された最大の差は 1.2e-13 で、shape がごく小さい (例: 0.001) ガンマ乱数で生じます。そこでは exp(log(u)/shape) が誤差を 1/shape 倍に拡大します。一様乱数の個数は、各ケースの後に runif(1) を引いて確かめます。採択・棄却の判定が 1 回でも異なれば、そこに現れます。
R の既定以外の RNG の種類 (rngRounding、rngLecuyer、rngStreams) も同じ方法で確かめ、さらに各ケースの後の状態 (getRandomSeed() と R の .Random.seed。種類のコードと、L’Ecuyer-CMRG では 7 個の整数すべて) を比べます。L’Ecuyer-CMRG の一様乱数は double による厳密な整数演算で計算するので、runif と離散値はすべて R と完全に一致します。連続値 (rnorm、rexp、rgamma) は上と同じく 1e-12 以内で比べます。
method、統計量の名前など) は完全に一致しなければなりません。R がエラーを出す入力では、エラーを投げなければなりません。| 桁落ちで決まるフィールドは、スケールの下限を設けて | 差 | / max( | 期待値 | , scale) で測ります。t 検定の信頼区間では | mu | を、Wilcoxon の信頼区間と推定値では、データの絶対値の最大と | mu | (真の根が 0 のとき) を使います。 |
ansariTest の正規近似の限界と推定値は、R でも rStats.js でも許容誤差 1e-4 の uniroot で求めた根なので、許容誤差 1e-4 で比較します。同じ根の探索を同じ階段関数に対して行うので、R とは 3e-15 以内で一致します。R が “object ‘cint’ not found” で失敗するケース (24 件、differences-from-r.md を参照) は、oracle が exact = FALSE で繰り返したケースと比較します。chisqTestSim、fisherTestSim) は完全に一致しなければなりません。各ケースは検定の直前に setSeed(seed) を呼び、その後に引く runif(1) が R と等しいことで、消費した一様乱数の個数を確かめます。chisqTest(x, y) は水準をコードポイント順に並べますが、R はロケールの照合順序を使いました。順序が異なるフィクスチャの 1 ケースでは、R の行と列を並べ替えてから比較します。| 行列とベクトルはノルムで比較します ( | 差 | / そのフィールドの絶対値の最大)。0 に近い要素は桁落ちで決まるため、要素ごとには比較できません。crossprod は列ノルムの積をスケールに使います。 |
args.cond に保存します。許容する誤差は scale × max(tolerance, 10 ε κ^k) で、k は既知の感度から決めます。chol は 1/2、solve、行列式、rcond、kappa、qr、ginv は 1、最小二乗の係数 (qrCoef、qrSolve) は 2 です。最小二乗の係数は、残差が 0 でないとき誤差が κ² に比例するためです。gap) で許容誤差を広げます。重複した値 (相対間隔が 1e-8 未満) と 0 の特異値に対するベクトルは一意に定まらないので、比較しません。qr: 数値ランク (最大の特異値の 1e-13 倍を超える特異値の個数) より後ろの列の Householder ベクトルは丸めだけで決まるため、比較しません。LAPACK モードでランク落ちした行列では、Q、qty、qy、係数も比較しません。test/linalg-invariants.test.ts で確かめます。qr と chol のピボットとランクは、R の限定的な列ピボット選択も含めて R と完全に一致します。info ブロックに記録しています。テストしている範囲では、既知の問題はありません。すべてのケースが、許容誤差の範囲で R または高精度の参照値と一致します。
その範囲の外では、次数が数千で x < ν の besselJ と besselY の相対誤差が最大で 1.5e-14 ほどになります (たとえば besselY(3000, 3949) で 1.04e-14、besselJ(1000, 1500) で 9.9e-15、besselJ(2000, 2999) で 1.5e-14。0.16.0 で、次数 1100〜9000 の 25 点を Rmpfr と比べて測りました)。漸化式と連分数の多数のステップで丸め誤差が積み重なるためです。
R の qnbinom は 13 個の入力で終了しません (例: qnbinom(-1, 0.5, 1e-10, log.p = TRUE))。oracle/nbinom.R は各ケースを 2 秒に制限し、それらの行をフィクスチャから外します。
ほとんどのカーネルは SciPy の xsf (BSD-3) と Boost.Math (BSL-1.0) から移植しています。線形代数は LAPACK のリファレンス実装 (修正 BSD)、LINPACK、EISPACK (パブリックドメイン) に従っています。THIRD_PARTY_NOTICES.md を参照してください。R に並ぶ、または R を上回るために加えた主な変更は次のとおりです。
src/special/loader.ts) と、指数分布、Weibull 分布、ロジスティック分布、対数正規分布における exp(−t) の引数 (src/internal/ddexp.ts) は double-double で組み立てます。t = 700 では t の丸めだけで 6e-14 の誤差になります。src/special/incbet.ts)。パラメータが大きく x が平均に近いとき、連分数は x の丸めにきわめて敏感です (a = 1e5, x = 0.999 で 6e-13)。係数の a + b と b − 1 も double-double にします。a + b を丸めると小さい b の桁が失われるためです (a = 1e8, b = 0.01 で 5e-11)。2^52 以上では 1 の刻みが失われるので、すべての係数を double-double にします。src/special/igam.ts)。a < 4 かつ 1.1 < x < 64 では、上側の裾の連分数の丸め誤差が double で約 120 eps まで積み重なっていました (pchisq(2.3569217160548437, 1, lower.tail = FALSE) が 45 eps ずれていた。R 4.6.2 はそこで 25 eps、その範囲で最大約 60 eps ずれる)。その範囲では漸化式を double-double で計算し (1 eps 以内)、pgamma と pchisq (df < 8) の上側の裾は数 eps 以内になりました (1.2.0)。1 回の呼び出しが最大で約 1 µs から約 8 µs になります。Temme の展開をその領域全体で double-double の逸脱度から (src/special/igam.ts)。xsf の領域 (20 < a < 200 で |
x − a | / a < 0.3、それより大きい a で 4.5 標準偏差以内) では、丸めた η から erfc と e^(−a η²/2) を作っていたため最大約 35 eps ずれていました (pgamma(256.71, 199, lower.tail = FALSE) が 22 eps、R は 0.2 eps)。小さいほうの裾を、その外側と同じく double-double の逸脱度から作るようにし、約 4 eps 以内になりました (1.2.0)。 |
src/special/loader.ts)。a < 0.01 では x^a e^−x / Γ(a + 1) を e^−x e^(a log x − lgamma1p(a)) で計算します。Loader の形では stirlerr(a) と log(2πa) が打ち消し合い、最大約 50 eps ずれていました (pgamma(100, 1e-100, lower.tail = FALSE) が 51 eps、R は 52 eps) (1.2.0)。src/special/incbetTemme.ts)。ab/(a + b) >= 25 で、x が展開の収束する範囲にあるとき、不完全ベータ関数は Temme の展開を使います。指数は Loader の二項密度と同じ double-double の逸脱度です。裾では、erfc の項と補正項の桁落ちを避けるために式を組み替えた形を使います (performance.md を参照)。src/special/incbetLargeA.ts)。a >= 15、b <= 40 で x が中央部にあるとき、u = −(a + (b − 1)/2) log x の不完全ガンマ関数による DiDonato と Morris の展開が連分数に代わります。前因子は Loader の二項密度と同じ double-double の逸脱度から作り、u は double-double で作り、先頭の不完全ガンマ関数は後ろから評価して double-double で仕上げる連分数で求めます (igam.ts の前向きの評価は小さな u では早く止まりすぎます)。そのため、これを通るケースは 1.2e-15 以内です (performance.md を参照)。src/special/lgammaDiff.ts)。a ≪ b のときに使い、lbeta、前置因子、級数の補数で利用します。lgamma の零点 (1 と 2) の近くでは lgamma1p を使います (xsf の式は x = 1 + 1e-10 で相対誤差 6e-7 になります)。分位点のソルバー (src/internal/invert.ts)。u = log x で Newton/Halley 反復を行い、 |
log x | が大きいときは x そのもので仕上げます (そのときは u での刻みが x の精度を制限するため)。丸めのために動かない刻みは収束とみなします。qlnorm は、z の丸めが sdlog |
z | 倍に拡大されるため、正規分布の分位点に double-double の Newton 補正を 1 回かけます。qgeom は閉じた式を使います。離散分布の二分探索は、答えが 2^53 を超えたら止めます。 |
pnorm の裾では、exp(−q²/2) を正確に計算するために q = hi + lo (hi は 1/16 刻みの格子上) に分けます。logP が −2048 未満 (Cephes の裾の近似の範囲外) の qnorm は、pnorm を使った Newton ステップを 1 回加えます。非心分布 (src/internal/poissonMixture.ts)。ポアソン混合は最頻値から外側へ向かって正の項だけで合計し、漸化式は double-double で、指数は別に保持します (仮数 × 2^(500k))。そのため対数での答えは 1e-300 未満でも正確です。非心 t は、片側の二重指数型の変数変換を施したあとの正の被積分関数を積分するので、t と δ の符号が逆でも桁落ちが起きません。遠い裾では部分積分した形を使い、さらに先では |
t | ^(−ν) のべき乗則を使います。非心カイ二乗分布の密度は、x × ncp が非常に大きいときベッセル関数の漸近展開を使います。 |
src/internal/summary.ts)。平均、分散、相関は double-double で合計し、偏差は掛け算の前に 2 のべき乗で尺度を変えるので、1e±200 付近のデータでもオーバーフローもアンダーフローもしません。これにかかる速度面のコストと、double-double の効果を測る方法は performance.md で説明しています。