`logDensity` for the gamma distribution is off by about 19 orders of magnitude at `k = 1e300`, `θ = 1e-300`

未关闭
  • 哪里不对？
#236 0 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

还没有人认领这个 Issue。

评估

难度
4/5
预计耗时
3-5 天
新手友好度
48/100
Issue 类型
缺陷
描述清晰度
基本清楚
活跃度
活跃
技术栈
haskell
领域
backend

调研方向

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
30 天内没有已合并 PR

环境准备

这个项目没有提供开发容器、Dockerfile 或贡献指南，环境需要你自己搭建：先看它的 README，通用步骤见我们的新手贡献指南。

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

haskell/statistics 的其他 Issue

查看 haskell/statistics 的全部 Issue

相似的 Issue

更多 Haskell Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。