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

Đang mở
  • Có gì không ổn?
#236 0 bình luận 0 reaction 0 người được giao Xem trên GitHub

Chưa có ai nhận issue này.

Đánh giá

Độ khó
4/5
Thời gian dự kiến
3-5 ngày
Mức phù hợp với người mới
48/100
Loại issue
Lỗi
Độ rõ ràng
Khá rõ ràng
Mức độ hoạt động
Sôi nổi
Công nghệ
haskell
Lĩnh vực
backend

Hướng nghiên cứu

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.

Do mô hình lập chỉ mục viết ra từ nội dung của issue.

Mô tả

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!

Ngôn ngữ chính
Haskell
Star
309
Fork
77
Không có pull request nào được merge trong 30 ngày

Chuẩn bị môi trường

Dự án này không cung cấp dev container, Dockerfile hay hướng dẫn đóng góp, nên bạn cần tự thiết lập môi trường: hãy bắt đầu từ README và xem hướng dẫn đóng góp lần đầu của chúng tôi để biết các bước chung.

Bắt đầu từ đâu

  1. Đọc hết issue, rồi đọc hướng dẫn đóng góp của dự án.
  2. Bình luận trên issue rằng bạn sẽ nhận — tránh hai người làm cùng một việc.
  3. Fork repository và làm thay đổi trên một nhánh.
  4. Mở pull request có tham chiếu số hiệu của issue.

Issue khác của haskell/statistics

Tất cả issue của haskell/statistics

Issue tương tự

Thêm issue về Haskell

Nhận issue mới trong hộp thư của bạn

Bản tóm tắt ngắn những issue GitHub phù hợp với người mới.