`logDensity` for the gamma distribution is off by about 19 orders of magnitude at `k = 1e300`, `θ = 1e-300`
まだ誰も着手していません。
評価
調査の方向性
Start by reading Statistics.Distribution.Gamma.logDensity and the kbn summation it uses; the issue provides a minimal Main.hs reproduction runnable with runghc. Investigate how to combine the large terms accurately for the stated shape and scale, then verify that the result for gammaDistr 1e300 1e-300 at 1 matches the reported reference value.
索引モデルが issue の本文から書いたものです。
説明
Hi! First of all, thanks a lot for your awesome library!
Unfortunately, as said in the title: logDensity of the gamma distribution with shape 1e300 and scale 1e-300 (both positive, so inside the documented domain), evaluated at 1, returns -6.111639493843688e286. The correct value at these doubles is -3.0080755755517465e+267. Both are negative, but they are about 19 orders of magnitude apart.
To reproduce it, save this as Main.hs and run runghc Main.hs with statistics-0.16.5.0 exposed (tested with math-functions-0.3.4.4, on GHC 9.6.7 and 9.14.1):
import Statistics.Distribution (logDensity)
import Statistics.Distribution.Gamma (gammaDistr)
main :: IO ()
main = print (logDensity (gammaDistr 1e300 1e-300) 1)
which prints:
-6.111639493843688e286
Here is the computation of the correct value from the gamma log density, (k-1)*log x - x/θ - k*log θ - logGamma k, at the exact binary values of the doubles passed for k and θ (with mp.dps raised until the 17 printed digits stopped changing):
from mpmath import mp, mpf, log, loggamma
mp.dps = 200
k = mpf(1e300) # double nearest 1e300
theta = mpf(1e-300) # double nearest 1e-300
x = mpf(1.0)
v = (k-1)*log(x) - x/theta - k*log(theta) - loggamma(k)
# v = -3.0080755755517465e+267
The double nearest 1e300 times the double nearest 1e-300 is 1 + 7.7564e-17, so x/θ and k differ by about 7.8e283 (at the exact 10^300 and 10^-300 they would be equal, and the formula gives +344.47, which no Double input can produce).
The terms log x * (k-1), -(x/θ), -logGamma k and -log θ * k are summed in Gamma.logDensity with kbn, the Kahan-Babuška-Neumaier sum from math-functions, but for these parameters they have magnitudes up to about 7e302, so even kbn can't recover the exact result in Double. Would it be possible to fix this? (combining the -x/θ and -k*log θ terms with the asymptotic expansion of logGamma analytically before evaluating would avoid summing terms about 2e35 times larger than the result)
Thanks in advance!
- 主要言語
- Haskell
- スター
- 309
- フォーク
- 77
- PR マージ指標
- 30日以内にマージされた PR はありません
環境構築
このプロジェクトには開発コンテナ、Dockerfile、コントリビューションガイドがありません。まず README を読み、一般的な手順ははじめてのコントリビューションガイドを参照してください。
はじめの一歩
- issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
- 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
- リポジトリをフォークし、ブランチを切って変更します。
- issue 番号を参照したプルリクエストを送ります。
haskell/statistics のほかの issue
-
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
haskell/statistics#235 ·
-
難易度 4/5 3〜5日 初心者へのやさしさ 48/100
haskell/statistics#237 ·
-
難易度 3/5 1〜2日 初心者へのやさしさ 58/100
haskell/statistics#234 ·
-
難易度 4/5 3〜5日 初心者へのやさしさ 42/100
haskell/statistics#231 · コメント 3 件 ·
-
難易度 4/5 3〜5日 初心者へのやさしさ 30/100
haskell/statistics#215 ·
haskell/statistics の issue をすべて見る
似ている issue
-
難易度 2/5 1〜3時間 初心者へのやさしさ 72/100
メンテナーはふだん 1 日以内に返信
-
Bug
難易度 2/5 1〜3時間 初心者へのやさしさ 68/100
-
pkg:crucible-mir type:bug
難易度 2/5 1〜3時間 初心者へのやさしさ 75/100
-
Heroes in statisticsオープン
難易度 2/5 1〜3時間 初心者へのやさしさ 68/100
-
bug
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
GaloisInc/llvm-pretty#211 ·