`logDensity` returns `0.0` for a beta distribution with shape parameters `1e200` at `0.5`
还没有人认领这个 Issue。
评估
调研方向
Start by locating logDensity for the beta distribution and the logBeta implementation used by math-functions; the report identifies cancellation between these calculations for shape parameters 1e200 at 0.5. Check existing beta distribution tests and add coverage for this reproduction. Done means the computed value is close to 230.37929153703981 rather than 0.0.
由索引模型根据 Issue 内容生成。
描述
Hello,
Unfortunately, as said in the title: logDensity returns 0.0 for a beta distribution with both shape parameters equal to 1e200, evaluated at 0.5. To reproduce it, save the following as Main.hs and run runghc Main.hs with statistics-0.16.5.0 exposed:
import qualified Statistics.Distribution as D
import Statistics.Distribution.Beta (betaDistr)
main :: IO ()
main = print (D.logDensity (betaDistr 1e200 1e200) 0.5)
It prints:
0.0
I get this on both GHC 9.6.7 and 9.14.1, with statistics-0.16.5.0 and math-functions-0.3.4.4.
The value should be approximately 230.37929153703981. Here is how to check it with mpmath (700 decimal digits of precision, starting from the exact binary floating-point values passed for the shape parameters and x):
from mpmath import mp
mp.dps = 700
a = mp.mpf(1e200)
x = mp.mpf(0.5)
log_beta = 2 * mp.loggamma(a) - mp.loggamma(2 * a)
log_density = 2 * (a - 1) * mp.log(x) - log_beta
print(mp.nstr(log_density, 17)) # 230.37929153703981
Both shape parameters are positive, so the input is inside the documented domain. The problem comes from the formula (a-1)*log x + (b-1)*log1p(-x) - logBeta a b used in logDensity. Here a = b = 1e200 and x = 0.5. So the (a-1) and (b-1) terms sum to 2*(1e200-1)*log 0.5 = -1.3862943611e200, and logBeta 1e200 1e200 (from the allStirling path in math-functions) is -1.3862943611e200 too. Both are negative, but their true difference is about 230 (the log density). In Double they cancel out, so the subtraction is 0.0 and can't recover a single digit. Would it be possible to fix this? (combining the (a-1)*log x and logBeta terms analytically before evaluating would do it, so they don't cancel out)
Thanks in advance!
- 主要语言
- Haskell
- 星标
- 309
- 派生
- 77
- PR 合并指标
- 30 天内没有已合并 PR
环境准备
这个项目没有提供开发容器、Dockerfile 或贡献指南，环境需要你自己搭建：先看它的 README，通用步骤见我们的新手贡献指南。
从这里开始
- 先读完整个 Issue，再读项目的贡献指南。
- 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
- Fork 仓库，在一个分支上完成修改。
- 提交 Pull Request，并在描述里引用这个 Issue 编号。
haskell/statistics 的其他 Issue
-
难度 2/5 1-3 小时 新手友好度 78/100
haskell/statistics#235 ·
-
难度 4/5 3-5 天 新手友好度 48/100
haskell/statistics#236 ·
-
难度 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
-
难度 2/5 1-3 小时 新手友好度 68/100
-
bug
难度 2/5 1-3 小时 新手友好度 78/100
GaloisInc/llvm-pretty#211 ·