-- |
-- 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 :: Double -> Double -> Double
probability Double
0      Double
0     = Double
1
probability Double
0      Double
1     = Double
0
probability Double
lambda Double
x
  | forall a. RealFloat a => a -> Bool
isInfinite Double
lambda    = Double
0
  | Double
x forall a. Ord a => a -> a -> Bool
< Double
0                = Double
0
  | Double
x forall a. Ord a => a -> a -> Bool
<= Double
lambda forall a. Num a => a -> a -> a
* Double
m_tiny = forall a. Floating a => a -> a
exp (-Double
lambda)
  | Double
lambda forall a. Ord a => a -> a -> Bool
< Double
x forall a. Num a => a -> a -> a
* Double
m_tiny  = forall a. Floating a => a -> a
exp (-Double
lambda forall a. Num a => a -> a -> a
+ Double
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
lambda forall a. Num a => a -> a -> a
- Double -> Double
logGamma (Double
xforall a. Num a => a -> a -> a
+Double
1))
  | Bool
otherwise            = forall a. Floating a => a -> a
exp (-(Double -> Double
stirlingError Double
x) forall a. Num a => a -> a -> a
- Double -> Double -> Double
bd0 Double
x Double
lambda) forall a. Fractional a => a -> a -> a
/
                           (Double
m_sqrt_2_pi forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
x)

-- -- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy
-- -- of the Poisson Law".  This function is unused because 'directEntropy'
-- -- 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 :: Double -> [Double]
powers Double
x = forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (\Double
y -> forall a. a -> Maybe a
Just (Double
yforall a. Num a => a -> a -> a
*Double
x,Double
yforall a. Num a => a -> a -> a
*Double
x)) Double
1

-- | Returns an upper bound according to theorem 2 of "Sharp Bounds on
-- the Entropy of the Poisson Law"
alyThm2Upper :: Double -> [Double] -> Double
alyThm2Upper :: Double -> [Double] -> Double
alyThm2Upper Double
lambda [Double]
coefficients =
  Double
1.4189385332046727 forall a. Num a => a -> a -> a
+ Double
0.5 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
lambda forall a. Num a => a -> a -> a
+
  Double -> [Double] -> Double
zipCoefficients Double
lambda [Double]
coefficients

-- | Returns the average of the upper and lower bounds according to
-- theorem 2.
alyThm2 :: Double -> [Double] -> [Double] -> Double
alyThm2 :: Double -> [Double] -> [Double] -> Double
alyThm2 Double
lambda [Double]
upper [Double]
lower =
  Double -> [Double] -> Double
alyThm2Upper Double
lambda [Double]
upper forall a. Num a => a -> a -> a
+ Double
0.5 forall a. Num a => a -> a -> a
* (Double -> [Double] -> Double
zipCoefficients Double
lambda [Double]
lower)

zipCoefficients :: Double -> [Double] -> Double
zipCoefficients :: Double -> [Double] -> Double
zipCoefficients Double
lambda [Double]
coefficients =
  (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a. Num a => a -> a -> a
(*)) (forall a b. [a] -> [b] -> [(a, b)]
zip (Double -> [Double]
powers forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => a -> a
recip Double
lambda) [Double]
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 :: [Double]
upperCoefficients4 = [Double
1forall a. Fractional a => a -> a -> a
/Double
12, Double
1forall a. Fractional a => a -> a -> a
/Double
24, -Double
103forall a. Fractional a => a -> a -> a
/Double
180, -Double
13forall a. Fractional a => a -> a -> a
/Double
40, -Double
1forall a. Fractional a => a -> a -> a
/Double
210]

lowerCoefficients4 :: [Double]
lowerCoefficients4 :: [Double]
lowerCoefficients4 = [Double
0,Double
0,Double
0, -Double
105forall a. Fractional a => a -> a -> a
/Double
4, -Double
210, -Double
2275forall a. Fractional a => a -> a -> a
/Double
18, -Double
167forall a. Fractional a => a -> a -> a
/Double
21, -Double
1forall a. Fractional a => a -> a -> a
/Double
72]

upperCoefficients6 :: [Double]
upperCoefficients6 :: [Double]
upperCoefficients6 = [Double
1forall a. Fractional a => a -> a -> a
/Double
12, Double
1forall a. Fractional a => a -> a -> a
/Double
24, Double
19forall a. Fractional a => a -> a -> a
/Double
360, Double
9forall a. Fractional a => a -> a -> a
/Double
80, -Double
38827forall a. Fractional a => a -> a -> a
/Double
2520,
                      -Double
74855forall a. Fractional a => a -> a -> a
/Double
1008, -Double
73061forall a. Fractional a => a -> a -> a
/Double
2520, -Double
827forall a. Fractional a => a -> a -> a
/Double
720, -Double
1forall a. Fractional a => a -> a -> a
/Double
990]

lowerCoefficients6 :: [Double]
lowerCoefficients6 :: [Double]
lowerCoefficients6 = [Double
0,Double
0,Double
0,Double
0,Double
0, -Double
3465forall a. Fractional a => a -> a -> a
/Double
2, -Double
45045, -Double
466235forall a. Fractional a => a -> a -> a
/Double
4, -Double
531916forall a. Fractional a => a -> a -> a
/Double
9,
                      -Double
56287forall a. Fractional a => a -> a -> a
/Double
10, -Double
629forall a. Fractional a => a -> a -> a
/Double
11, -Double
1forall a. Fractional a => a -> a -> a
/Double
156]

upperCoefficients8 :: [Double]
upperCoefficients8 :: [Double]
upperCoefficients8 = [Double
1forall a. Fractional a => a -> a -> a
/Double
12, Double
1forall a. Fractional a => a -> a -> a
/Double
24, Double
19forall a. Fractional a => a -> a -> a
/Double
360, Double
9forall a. Fractional a => a -> a -> a
/Double
80, Double
863forall a. Fractional a => a -> a -> a
/Double
2520, Double
1375forall a. Fractional a => a -> a -> a
/Double
1008,
                      -Double
3023561forall a. Fractional a => a -> a -> a
/Double
2520, -Double
15174047forall a. Fractional a => a -> a -> a
/Double
720, -Double
231835511forall a. Fractional a => a -> a -> a
/Double
5940,
                      -Double
18927611forall a. Fractional a => a -> a -> a
/Double
1320, -Double
58315591forall a. Fractional a => a -> a -> a
/Double
60060, -Double
23641forall a. Fractional a => a -> a -> a
/Double
3640,
                      -Double
1forall a. Fractional a => a -> a -> a
/Double
2730]

lowerCoefficients8 :: [Double]
lowerCoefficients8 :: [Double]
lowerCoefficients8 = [Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0, -Double
2027025forall a. Fractional a => a -> a -> a
/Double
8, -Double
15315300, -Double
105252147,
                      -Double
178127950, -Double
343908565forall a. Fractional a => a -> a -> a
/Double
4, -Double
10929270, -Double
3721149forall a. Fractional a => a -> a -> a
/Double
14,
                      -Double
7709forall a. Fractional a => a -> a -> a
/Double
15, -Double
1forall a. Fractional a => a -> a -> a
/Double
272]

upperCoefficients10 :: [Double]
upperCoefficients10 :: [Double]
upperCoefficients10 = [Double
1forall a. Fractional a => a -> a -> a
/Double
12, Double
1forall a. Fractional a => a -> a -> a
/Double
24, Double
19forall a. Fractional a => a -> a -> a
/Double
360, Double
9,Double
80, Double
863forall a. Fractional a => a -> a -> a
/Double
2520, Double
1375forall a. Fractional a => a -> a -> a
/Double
1008,
                       Double
33953forall a. Fractional a => a -> a -> a
/Double
5040, Double
57281forall a. Fractional a => a -> a -> a
/Double
1440, -Double
2271071617forall a. Fractional a => a -> a -> a
/Double
11880,
                       -Double
1483674219forall a. Fractional a => a -> a -> a
/Double
176, -Double
31714406276557forall a. Fractional a => a -> a -> a
/Double
720720,
                       -Double
7531072742237forall a. Fractional a => a -> a -> a
/Double
131040, -Double
1405507544003forall a. Fractional a => a -> a -> a
/Double
65520,
                       -Double
21001919627forall a. Fractional a => a -> a -> a
/Double
10080, -Double
1365808297forall a. Fractional a => a -> a -> a
/Double
36720,
                       -Double
26059forall a. Fractional a => a -> a -> a
/Double
544, -Double
1forall a. Fractional a => a -> a -> a
/Double
5814]

lowerCoefficients10 :: [Double]
lowerCoefficients10 :: [Double]
lowerCoefficients10 = [Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,-Double
130945815forall a. Fractional a => a -> a -> a
/Double
2, -Double
7638505875,
                       -Double
438256243425forall a. Fractional a => a -> a -> a
/Double
4, -Double
435477637540, -Double
3552526473925forall a. Fractional a => a -> a -> a
/Double
6,
                       -Double
857611717105forall a. Fractional a => a -> a -> a
/Double
3, -Double
545654955967forall a. Fractional a => a -> a -> a
/Double
12, -Double
5794690528forall a. Fractional a => a -> a -> a
/Double
3,
                       -Double
578334559forall a. Fractional a => a -> a -> a
/Double
42, -Double
699043forall a. Fractional a => a -> a -> a
/Double
133, -Double
1forall a. Fractional a => a -> a -> a
/Double
420]

upperCoefficients12 :: [Double]
upperCoefficients12 :: [Double]
upperCoefficients12 = [Double
1forall a. Fractional a => a -> a -> a
/Double
12, Double
1forall a. Fractional a => a -> a -> a
/Double
24, Double
19forall a. Fractional a => a -> a -> a
/Double
360, Double
863forall a. Fractional a => a -> a -> a
/Double
2520, Double
1375forall a. Fractional a => a -> a -> a
/Double
1008,
                       Double
33953forall a. Fractional a => a -> a -> a
/Double
5040, Double
57281forall a. Fractional a => a -> a -> a
/Double
1440, Double
3250433forall a. Fractional a => a -> a -> a
/Double
11880,
                       Double
378351forall a. Fractional a => a -> a -> a
/Double
176, -Double
37521922090657forall a. Fractional a => a -> a -> a
/Double
720720,
                       -Double
612415657466657forall a. Fractional a => a -> a -> a
/Double
131040, -Double
3476857538815223forall a. Fractional a => a -> a -> a
/Double
65520,
                       -Double
243882174660761forall a. Fractional a => a -> a -> a
/Double
1440, -Double
34160796727900637forall a. Fractional a => a -> a -> a
/Double
183600,
                       -Double
39453820646687forall a. Fractional a => a -> a -> a
/Double
544, -Double
750984629069237forall a. Fractional a => a -> a -> a
/Double
81396,
                       -Double
2934056300989forall a. Fractional a => a -> a -> a
/Double
9576, -Double
20394527513forall a. Fractional a => a -> a -> a
/Double
12540,
                       -Double
3829559forall a. Fractional a => a -> a -> a
/Double
9240, -Double
1forall a. Fractional a => a -> a -> a
/Double
10626]

lowerCoefficients12 :: [Double]
lowerCoefficients12 :: [Double]
lowerCoefficients12 = [Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,Double
0,
                       -Double
105411381075forall a. Fractional a => a -> a -> a
/Double
4, -Double
5270569053750, -Double
272908057767345forall a. Fractional a => a -> a -> a
/Double
2,
                       -Double
1051953238104769, -Double
24557168490009155forall a. Fractional a => a -> a -> a
/Double
8,
                       -Double
3683261873403112, -Double
5461918738302026forall a. Fractional a => a -> a -> a
/Double
3,
                       -Double
347362037754732, -Double
2205885452434521forall a. Fractional a => a -> a -> a
/Double
100,
                       -Double
12237195698286forall a. Fractional a => a -> a -> a
/Double
35, -Double
16926981721forall a. Fractional a => a -> a -> a
/Double
22,
                       -Double
6710881forall a. Fractional a => a -> a -> a
/Double
155, -Double
1forall a. Fractional a => a -> a -> a
/Double
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 :: Double -> Double
directEntropy Double
lambda =
  forall a. Num a => a -> a
negate forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$
  forall a. (a -> Bool) -> [a] -> [a]
takeWhile (forall a. Ord a => a -> a -> Bool
< forall a. Num a => a -> a
negate Double
m_epsilon forall a. Num a => a -> a -> a
* Double
lambda) forall a b. (a -> b) -> a -> b
$
  forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall a. Ord a => a -> a -> Bool
< forall a. Num a => a -> a
negate Double
m_epsilon forall a. Num a => a -> a -> a
* Double
lambda)) forall a b. (a -> b) -> a -> b
$
  [ let x :: Double
x = Double -> Double -> Double
probability Double
lambda Double
k in Double
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
x | Double
k <- [Double
0..]]

-- | Compute the entropy of a Poisson distribution using the best available
-- method.
poissonEntropy :: Double -> Double
poissonEntropy :: Double -> Double
poissonEntropy Double
lambda
  | Double
lambda forall a. Eq a => a -> a -> Bool
== Double
0 = Double
0
  | Double
lambda forall a. Ord a => a -> a -> Bool
<= Double
10 = Double -> Double
directEntropy Double
lambda
  | Double
lambda forall a. Ord a => a -> a -> Bool
<= Double
12 = Double -> [Double] -> [Double] -> Double
alyThm2 Double
lambda [Double]
upperCoefficients4 [Double]
lowerCoefficients4
  | Double
lambda forall a. Ord a => a -> a -> Bool
<= Double
18 = Double -> [Double] -> [Double] -> Double
alyThm2 Double
lambda [Double]
upperCoefficients6 [Double]
lowerCoefficients6
  | Double
lambda forall a. Ord a => a -> a -> Bool
<= Double
24 = Double -> [Double] -> [Double] -> Double
alyThm2 Double
lambda [Double]
upperCoefficients8 [Double]
lowerCoefficients8
  | Double
lambda forall a. Ord a => a -> a -> Bool
<= Double
30 = Double -> [Double] -> [Double] -> Double
alyThm2 Double
lambda [Double]
upperCoefficients10 [Double]
lowerCoefficients10
  | Bool
otherwise = Double -> [Double] -> [Double] -> Double
alyThm2 Double
lambda [Double]
upperCoefficients12 [Double]
lowerCoefficients12