search

haskell-numerics / random-fu

A suite of Haskell libraries for representing, manipulating, and sampling random variables
42 stars 21 forks source link

Binomial CDF not numericaly stable #35

Closed rpglover64 closed 8 years ago

rpglover64 commented 8 years ago

integralBinomialCDF 5000 (1/52) 170 computes to Infinity, integralBinomialCDF 5000 (1/52) 200 computes to NaN.

idontgetoutmuch commented 8 years ago

Thanks for pointing this out. R has

> sum(dbinom(1:151,5000,1/52))
[1] 0.9999999

> sum(dbinom(1:152,5000,1/52))
[1] 1

I think the solution is to use

*Main> Data.List.sum $ map exp $ map (integralBinomialLogPdf 5000 (1/52)) (map fromIntegral [0..170])
0.9999999999978313
*Main> Data.List.sum $ map exp $ map (integralBinomialLogPdf 5000 (1/52)) (map fromIntegral [0..200])
0.9999999999999999

Any views?

I should add that this uses the method given in Catherine Loader's paper, the only reference for which seems to be in here:

https://github.com/scipy/scipy/issues/1147

idontgetoutmuch commented 8 years ago

I just made the change.