Skip to content
Browse files

Share the fast accurate Poisson between Gamma and Poisson.

  • Loading branch information...
1 parent 06a482c commit 9cd9058379bf91ffd3ccd54163f963cd4d1e8f34 @bos committed Aug 25, 2011
View
5 Statistics/Constants.hs
@@ -13,6 +13,7 @@ module Statistics.Constants
(
m_epsilon
, m_huge
+ , m_tiny
, m_1_sqrt_2
, m_2_sqrt_pi
, m_ln_sqrt_2_pi
@@ -29,6 +30,10 @@ m_huge :: Double
m_huge = 1.7976931348623157e308
{-# INLINE m_huge #-}
+m_tiny :: Double
+m_tiny = 2.2250738585072014e-308
+{-# INLINE m_tiny #-}
+
-- | The largest 'Int' /x/ such that 2**(/x/-1) is approximately
-- representable as a 'Double'.
m_max_exp :: Int
View
7 Statistics/Distribution/Gamma.hs
@@ -26,7 +26,8 @@ module Statistics.Distribution.Gamma
import Data.Typeable (Typeable)
import Statistics.Constants (m_huge, m_pos_inf, m_NaN)
-import Statistics.Math (incompleteGamma, pois)
+import Statistics.Distribution.Poisson.Internal as Poisson
+import Statistics.Math (incompleteGamma)
import qualified Statistics.Distribution as D
-- | The gamma distribution.
@@ -68,8 +69,8 @@ density (GD a l) x
| x <= 0 = 0
| a == 0 = if x == 0 then m_pos_inf else 0
| x == 0 = if a < 1 then m_pos_inf else if a > 1 then 0 else 1/l
- | a < 1 = (pois (x/l) a)*a/x
- | otherwise = (pois (x/l) (a-1))/l
+ | a < 1 = Poisson.probability (x/l) a * a / x
+ | otherwise = Poisson.probability (x/l) (a-1) / l
{-# INLINE density #-}
cumulative :: GammaDistribution -> Double -> Double
View
25 Statistics/Distribution/Poisson.hs
@@ -20,11 +20,13 @@ module Statistics.Distribution.Poisson
, poisson
-- * Accessors
, poissonLambda
+ -- * References
+ -- $references
) where
import Data.Typeable (Typeable)
import qualified Statistics.Distribution as D
-import Statistics.Math (logGamma, factorial)
+import qualified Statistics.Distribution.Poisson.Internal as I
newtype PoissonDistribution = PD {
poissonLambda :: Double
@@ -35,7 +37,8 @@ instance D.Distribution PoissonDistribution where
{-# INLINE cumulative #-}
instance D.DiscreteDistr PoissonDistribution where
- probability = probability
+ probability (PD lambda) x = I.probability lambda (fromIntegral x)
+ {-# INLINE probability #-}
instance D.Variance PoissonDistribution where
variance = poissonLambda
@@ -45,19 +48,15 @@ instance D.Mean PoissonDistribution where
mean = poissonLambda
{-# INLINE mean #-}
--- | Create poisson distribution.
+-- | Create Poisson distribution.
poisson :: Double -> PoissonDistribution
poisson l
- | l <= 0 =
- error $ "Statistics.Distribution.Poisson.poisson: lambda must be positive. Got " ++ show l
+ | l <= 0 = error $ "Statistics.Distribution.Poisson.poisson:\
+ \ lambda must be positive. Got " ++ show l
| otherwise = PD l
{-# INLINE poisson #-}
-probability :: PoissonDistribution -> Int -> Double
-probability (PD l) n
- | n < 0 = 0
- | l < 20 && n <= 100 = exp (-l) * l ** x / factorial n
- | otherwise = exp (x * log l - logGamma (x + 1) - l)
- where
- x = fromIntegral n
-{-# INLINE probability #-}
+-- $references
+--
+-- * Loader, C. (2000) Fast and Accurate Computation of Binomial
+-- Probabilities. <http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf>
View
32 Statistics/Distribution/Poisson/Internal.hs
@@ -0,0 +1,32 @@
+-- |
+-- 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
+ ) where
+
+import Statistics.Constants (m_sqrt_2_pi, m_tiny)
+import Statistics.Math (bd0, logGamma, stirlingError)
+
+-- | 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)
+{-# INLINE probability #-}
View
18 Statistics/Math.hs
@@ -29,10 +29,11 @@ module Statistics.Math
, logGammaL
-- ** Logarithm
, log1p
+ -- ** Stirling's approximation
+ , stirlingError
+ , bd0
-- * References
-- $references
- -- ** Fast Poisson
- , pois
) where
import Data.Int (Int64)
@@ -405,19 +406,6 @@ bd0 x np
s' | s' == s -> s'
| otherwise -> loop (j+1) (ej*vv) s'
--- | fast accurate poisson distribution via Catherine Loader's algorithm.
-pois :: Double -> Double -> Double
-pois lambda x
- | lambda == 0 = if x == 0 then 1 else 0
- | isInfinite lambda = 0
- | x < 0 = 0
- | x <= lambda * dbl_min = exp (-lambda)
- | lambda < x * dbl_min = exp (-lambda + x*(log lambda) - (logGamma (x+1)))
- | otherwise = exp (-(stirlingError x)-(bd0 x lambda) ) / (m_sqrt_2_pi * (sqrt x))
- where
- dbl_min = 2 ** (- 1021) :: Double
-
-
-- $references
--
-- * Broucke, R. (1973) Algorithm 446: Ten subroutines for the
View
1 statistics.cabal
@@ -56,6 +56,7 @@ library
Statistics.Test.NonParametric
Statistics.Types
other-modules:
+ Statistics.Distribution.Poisson.Internal
Statistics.Internal
build-depends:
aeson,

0 comments on commit 9cd9058

Please sign in to comment.
Something went wrong with that request. Please try again.