1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
-- |
-- 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 '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 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 according 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
|
