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