module Statistics.Distribution.Poisson.Internal
(
probability, poissonEntropy
) where
import Data.List (unfoldr)
import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny, m_epsilon)
import Numeric.SpecFunctions (logGamma, stirlingError )
import Numeric.SpecFunctions.Extra (bd0)
probability :: Double -> Double -> Double
probability 0 0 = 1
probability 0 1 = 0
probability lambda x
| isInfinite lambda = 0
| x < 0 = 0
| x <= lambda * m_tiny = exp (lambda)
| lambda < x * m_tiny = exp (lambda + x * log lambda logGamma (x+1))
| otherwise = exp ((stirlingError x) bd0 x lambda) /
(m_sqrt_2_pi * sqrt x)
powers :: Double -> [Double]
powers x = unfoldr (\y -> Just (y*x,y*x)) 1
alyThm2Upper :: Double -> [Double] -> Double
alyThm2Upper lambda coefficients =
1.4189385332046727 + 0.5 * log lambda +
zipCoefficients lambda coefficients
alyThm2 :: Double -> [Double] -> [Double] -> Double
alyThm2 lambda upper lower =
alyThm2Upper lambda upper + 0.5 * (zipCoefficients lambda lower)
zipCoefficients :: Double -> [Double] -> Double
zipCoefficients lambda coefficients =
(sum $ map (uncurry (*)) (zip (powers $ recip lambda) coefficients))
upperCoefficients4 :: [Double]
upperCoefficients4 = [1/12, 1/24, 103/180, 13/40, 1/210]
lowerCoefficients4 :: [Double]
lowerCoefficients4 = [0,0,0, 105/4, 210, 2275/18, 167/21, 1/72]
upperCoefficients6 :: [Double]
upperCoefficients6 = [1/12, 1/24, 19/360, 9/80, 38827/2520,
74855/1008, 73061/2520, 827/720, 1/990]
lowerCoefficients6 :: [Double]
lowerCoefficients6 = [0,0,0,0,0, 3465/2, 45045, 466235/4, 531916/9,
56287/10, 629/11, 1/156]
upperCoefficients8 :: [Double]
upperCoefficients8 = [1/12, 1/24, 19/360, 9/80, 863/2520, 1375/1008,
3023561/2520, 15174047/720, 231835511/5940,
18927611/1320, 58315591/60060, 23641/3640,
1/2730]
lowerCoefficients8 :: [Double]
lowerCoefficients8 = [0,0,0,0,0,0,0, 2027025/8, 15315300, 105252147,
178127950, 343908565/4, 10929270, 3721149/14,
7709/15, 1/272]
upperCoefficients10 :: [Double]
upperCoefficients10 = [1/12, 1/24, 19/360, 9,80, 863/2520, 1375/1008,
33953/5040, 57281/1440, 2271071617/11880,
1483674219/176, 31714406276557/720720,
7531072742237/131040, 1405507544003/65520,
21001919627/10080, 1365808297/36720,
26059/544, 1/5814]
lowerCoefficients10 :: [Double]
lowerCoefficients10 = [0,0,0,0,0,0,0,0,0,130945815/2, 7638505875,
438256243425/4, 435477637540, 3552526473925/6,
857611717105/3, 545654955967/12, 5794690528/3,
578334559/42, 699043/133, 1/420]
upperCoefficients12 :: [Double]
upperCoefficients12 = [1/12, 1/24, 19/360, 863/2520, 1375/1008,
33953/5040, 57281/1440, 3250433/11880,
378351/176, 37521922090657/720720,
612415657466657/131040, 3476857538815223/65520,
243882174660761/1440, 34160796727900637/183600,
39453820646687/544, 750984629069237/81396,
2934056300989/9576, 20394527513/12540,
3829559/9240, 1/10626]
lowerCoefficients12 :: [Double]
lowerCoefficients12 = [0,0,0,0,0,0,0,0,0,0,0,
105411381075/4, 5270569053750, 272908057767345/2,
1051953238104769, 24557168490009155/8,
3683261873403112, 5461918738302026/3,
347362037754732, 2205885452434521/100,
12237195698286/35, 16926981721/22,
6710881/155, 1/600]
directEntropy :: Double -> Double
directEntropy lambda =
negate . sum $
takeWhile (< negate m_epsilon * lambda) $
dropWhile (not . (< negate m_epsilon * lambda)) $
[ let x = probability lambda k in x * log x | k <- [0..]]
poissonEntropy :: Double -> Double
poissonEntropy lambda
| lambda == 0 = 0
| lambda <= 10 = directEntropy lambda
| lambda <= 12 = alyThm2 lambda upperCoefficients4 lowerCoefficients4
| lambda <= 18 = alyThm2 lambda upperCoefficients6 lowerCoefficients6
| lambda <= 24 = alyThm2 lambda upperCoefficients8 lowerCoefficients8
| lambda <= 30 = alyThm2 lambda upperCoefficients10 lowerCoefficients10
| otherwise = alyThm2 lambda upperCoefficients12 lowerCoefficients12