-- | -- Module : Statistics.Distribution.Poisson.Internal -- Copyright : (c) 2011 Bryan O'Sullivan -- License : BSD3 -- -- Maintainer : bos@serpentine.com -- Stability : experimental -- Portability : portable -- -- Internal code for the Poisson distribution. 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, choose, logFactorial) import Numeric.SpecFunctions.Extra (bd0) -- | An unchecked, non-integer-valued version of Loader's saddle point -- algorithm. 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) -- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy -- of the Poisson Law". This function is unused because 'directEntorpy' -- is just as accurate and is faster by about a factor of 4. alyThm1 :: Double -> Double alyThm1 lambda = sum (takeWhile (\x -> abs x >= m_epsilon * lll) alySeries) + lll where lll = lambda * (1 - log lambda) alySeries = [ alyc k * exp (fromIntegral k * log lambda - logFactorial k) | k <- [2..] ] alyc :: Int -> Double alyc k = sum [ parity j * choose (k-1) j * log (fromIntegral j+1) | j <- [0..k-1] ] where parity j | even (k-j) = -1 | otherwise = 1 -- | Returns [x, x^2, x^3, x^4, ...] powers :: Double -> [Double] powers x = unfoldr (\y -> Just (y*x,y*x)) 1 -- | Returns an upper bound according to theorem 2 of "Sharp Bounds on -- the Entropy of the Poisson Law" alyThm2Upper :: Double -> [Double] -> Double alyThm2Upper lambda coefficients = 1.4189385332046727 + 0.5 * log lambda + zipCoefficients lambda coefficients -- | Returns the average of the upper and lower bounds accounding to -- theorem 2. 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)) -- Mathematica code deriving the coefficients below: -- -- poissonMoment[0, s_] := 1 -- poissonMoment[1, s_] := 0 -- poissonMoment[k_, s_] := -- Sum[s * Binomial[k - 1, j] * poissonMoment[j, s], {j, 0, k - 2}] -- -- upperSeries[m_] := -- Distribute[Integrate[ -- Sum[(-1)^(j - 1) * -- poissonMoment[j, \[Lambda]] / (j * (j - 1)* \[Lambda]^j), -- {j, 3, 2 m - 1}], -- \[Lambda]]] -- -- lowerSeries[m_] := -- Distribute[Integrate[ -- poissonMoment[ -- 2 m + 2, \[Lambda]] / ((2 m + -- 1)*\[Lambda]^(2 m + 2)), \[Lambda]]] -- -- upperBound[m_] := upperSeries[m] + (Log[2*Pi*\[Lambda]] + 1)/2 -- -- lowerBound[m_] := upperBound[m] + lowerSeries[m] 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] -- | Compute entropy directly from its definition. This is just as accurate -- as 'alyThm1' for lambda <= 1 and is faster, but is slow for large lambda, -- and produces some underestimation due to accumulation of floating point -- error. 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..]] -- | Compute the entropy of a poisson distribution using the best available -- method. 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