# Commits

committed b5b1eab

Improvements to Poisson entropy computation

• Participants
• Parent commits e09280e

# File Statistics/Distribution/Poisson/Internal.hs

` `
` import Data.List(unfoldr)`
` import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny, m_epsilon)`
`-import Numeric.SpecFunctions           (logGamma, stirlingError, choose)`
`+import Numeric.SpecFunctions (logGamma, stirlingError, choose, logFactorial)`
` import Numeric.SpecFunctions.Extra     (bd0)`
` `
` -- | An unchecked, non-integer-valued version of Loader's saddle point`
` {-# INLINE probability #-}`
` `
` -- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy`
`--- of the Poisson Law" by José A. Adell, Alberto Lekuona, and Yaming`
`--- Yu.  This is highly robust for lambda <= 1.`
`+-- 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 `
`-                          - logGamma (fromIntegral k))`
`+          [ alyc k * exp (fromIntegral k * log lambda - logFactorial k)`
`           | k <- [2..] ]`
` `
` alyc :: Int -> Double`
`           | even (k-j) = -1`
`           | otherwise  = 1`
`                          `
`-                         `
`--- | Compute entropy using Theorem 2 from the same paper. This is highly`
`--- accurate for large lambda.`
`-alyThm2 :: Double -> Double`
`-alyThm2 lambda =`
`+-- | 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 +`
`-  (sum \$ map (uncurry (*)) (zip powers coefficients))`
`-  where powers = unfoldr (\x -> Just (x/lambda,x/lambda)) 1`
`-        coefficients =`
`-          [1/12, 1/24, 19/360, 9/80, 863/2520, -123365/1008, -3023561/2520,`
`-           -808157/720, -984451/5940, -1151/440, -1/1716]`
`+  zipCoefficients lambda coefficients`
` `
`--- | Compute entropy by brute force. Limits on floating point precision`
`--- cause this to underestimate by quite a bit, but for intermediate values`
`--- of lambda it's the best we can do.`
`-bruteforce :: Double -> Double`
`-bruteforce lambda =   `
`+-- | 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. Use 'alyThm1' for`
`--- small lambda, 'alyThm2' for large lambda, and 'bruteforce' plus a`
`--- linear spline for intermediate lambda.`
`+-- | Compute the entropy of a poisson distribution using the best available`
`+-- method.`
` poissonEntropy :: Double -> Double`
` poissonEntropy lambda`
`   | lambda == 0 = 0`
`-  | lambda <= 1 = alyThm1 lambda`
`-  | lambda >= 4.671691395381868 = alyThm2 lambda`
`-  | otherwise = bruteforce lambda `
`-                - 7.314355647771469e-2 * lambda`
`-                + 0.3417041234245674`
`+  | 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`