statistics 0.10.4.1 → 0.10.5.0

raw patch · 20 files changed

+520/−107 lines, 20 filesdep ~mwc-randomnew-uploader

Dependency ranges changed: mwc-random

Files

• ChangeLog +22/−1
• Statistics/Distribution.hs +35/−1
• Statistics/Distribution/Beta.hs +14/−1
• Statistics/Distribution/Binomial.hs +20/−5
• Statistics/Distribution/CauchyLorentz.hs +6/−0
• Statistics/Distribution/ChiSquared.hs +13/−1
• Statistics/Distribution/Exponential.hs +14/−5
• Statistics/Distribution/FDistribution.hs +26/−7
• Statistics/Distribution/Gamma.hs +16/−2
• Statistics/Distribution/Geometric.hs +98/−18
• Statistics/Distribution/Hypergeometric.hs +14/−1
• Statistics/Distribution/Normal.hs +10/−4
• Statistics/Distribution/Poisson.hs +13/−2
• Statistics/Distribution/Poisson/Internal.hs +148/−3
• Statistics/Distribution/StudentT.hs +16/−5
• Statistics/Distribution/Transform.hs +10/−1
• Statistics/Distribution/Uniform.hs +6/−0
• Statistics/Sample.hs +1/−1
• statistics.cabal +2/−45
• tests/Tests/Distribution.hs +36/−4

@@ -1,4 +1,25 @@- Changes in 0.10.2.0+Changes in 0.10.5.0++  * Enthropy type class for distributions is added.++  * Probability and probability density of distribution is given in+    log domain too.++Changes in 0.10.4.1++  * Bugfix for GHC < 7.6++Changes in 0.10.4.0++  * Support for versions of GHC older than 7.2 is discontinued.++  * All datatypes now support 'Data.Binary' and 'GHC.Generics'.++Changes in 0.10.3.0++  * Bug fixes++Changes in 0.10.2.0    * Bugs in DCT and IDCT are fixed. 

@@ -21,6 +21,8 @@     , Mean(..)     , MaybeVariance(..)     , Variance(..)+    , MaybeEntropy(..)+    , Entropy(..)       -- ** Random number generation     , ContGen(..)     , DiscreteGen(..)@@ -65,20 +67,36 @@ class Distribution  d => DiscreteDistr d where     -- | Probability of n-th outcome.     probability :: d -> Int -> Double+    probability d = exp . logProbability d+    {-# INLINE probability #-} +    -- | Logarithm of probability of n-th outcome+    logProbability :: d -> Int -> Double+    logProbability d = log . probability d+    {-# INLINE logProbability #-} --- | Continuous probability distributuion++-- | Continuous probability distributuion.+--+--   Minimal complete definition is 'quantile' and either 'density' or+--   'logDensity'. class Distribution d => ContDistr d where     -- | Probability density function. Probability that random     -- variable /X/ lies in the infinitesimal interval     -- [/x/,/x+/δ/x/) equal to /density(x)/⋅δ/x/     density :: d -> Double -> Double+    density d = exp . logDensity d+    {-# INLINE density #-}      -- | Inverse of the cumulative distribution function. The value     -- /x/ for which P(/X/≤/x/) = /p/. If probability is outside     -- of [0,1] range function should call 'error'     quantile :: d -> Double -> Double +    -- | Natural logarithm of density.+    logDensity :: d -> Double -> Double+    logDensity d = log . density d+    {-# INLINE logDensity #-}   -- | Type class for distributions with mean. 'maybeMean' should return@@ -116,6 +134,22 @@     stdDev   :: d -> Double     stdDev = sqrt . variance +-- | Type class for distributions with entropy, meaning Shannon entropy+--   in the case of a discrete distribution, or differential entropy in the+--   case of a continuous one.  'maybeEntropy' should return 'Nothing' if+--   entropy is undefined for the chosen parameter values.+class (Distribution d) => MaybeEntropy d where+  -- | Returns the entropy of a distribution, in nats, if such is defined.+  maybeEntropy :: d -> Maybe Double+  +-- | Type class for distributions with entropy, meaning Shannon+--   entropy in the case of a discrete distribution, or differential+--   entropy in the case of a continuous one.  If the distribution has+--   well-defined entropy for all valid parameter values then it+--   should be an instance of this type class.+class (MaybeEntropy d) => Entropy d where+  -- | Returns the entropy of a distribution, in nats.+  entropy :: d -> Double  -- | Generate discrete random variates which have given --   distribution.

@@ -23,7 +23,8 @@ import Data.Binary (Binary) import Data.Data (Data, Typeable) import GHC.Generics (Generic)-import Numeric.SpecFunctions           (incompleteBeta, invIncompleteBeta, logBeta)+import Numeric.SpecFunctions (+  incompleteBeta, invIncompleteBeta, logBeta, digamma) import Numeric.MathFunctions.Constants (m_NaN) import qualified Statistics.Distribution as D @@ -81,6 +82,18 @@ instance D.MaybeVariance BetaDistribution where   maybeVariance = Just . D.variance   {-# INLINE maybeVariance #-}++instance D.Entropy BetaDistribution where+  entropy (BD a b) =+    logBeta a b +    - (a-1) * digamma a+    - (b-1) * digamma b+    + (a+b-2) * digamma (a+b)+  {-# INLINE entropy #-}+    +instance D.MaybeEntropy BetaDistribution where+  maybeEntropy = Just . D.entropy+  {-# INLINE maybeEntropy #-}  instance D.ContDistr BetaDistribution where   density (BD a b) x

@@ -27,7 +27,9 @@ import Data.Data (Data, Typeable) import GHC.Generics (Generic) import qualified Statistics.Distribution as D-import Numeric.SpecFunctions (choose)+import qualified Statistics.Distribution.Poisson.Internal as I+import Numeric.SpecFunctions (choose,incompleteBeta)+import Numeric.MathFunctions.Constants (m_epsilon)   -- | The binomial distribution.@@ -59,7 +61,15 @@     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance +instance D.Entropy BinomialDistribution where+  entropy (BD n p)+    | n == 0 = 0+    | n <= 100 = directEntropy (BD n p)+    | otherwise = I.poissonEntropy (fromIntegral n * p) +instance D.MaybeEntropy BinomialDistribution where+  maybeEntropy = Just . D.entropy+ -- This could be slow for big n probability :: BinomialDistribution -> Int -> Double probability (BD n p) k@@ -70,15 +80,13 @@  -- Summation from different sides required to reduce roundoff errors cumulative :: BinomialDistribution -> Double -> Double-cumulative d@(BD n _) x+cumulative (BD n p) x   | isNaN x      = error "Statistics.Distribution.Binomial.cumulative: NaN input"   | isInfinite x = if x > 0 then 1 else 0   | k <  0       = 0   | k >= n       = 1-  | k <  m       = D.sumProbabilities d 0 k-  | otherwise    = 1 - D.sumProbabilities d (k+1) n+  | otherwise    = incompleteBeta (fromIntegral (n-k)) (fromIntegral (k+1)) (1 - p)   where-    m = floor (mean d)     k = floor x {-# INLINE cumulative #-} @@ -89,6 +97,13 @@ variance :: BinomialDistribution -> Double variance (BD n p) = fromIntegral n * p * (1 - p) {-# INLINE variance #-}++directEntropy :: BinomialDistribution -> Double+directEntropy d@(BD n _) =   +  negate . sum $+  takeWhile (< negate m_epsilon) $+  dropWhile (not . (< negate m_epsilon)) $+  [ let x = probability d k in x * log x | k <- [0..n]]  -- | Construct binomial distribution. Number of trials must be --   non-negative and probability must be in [0,1] range

@@ -69,3 +69,9 @@  instance D.ContGen CauchyDistribution where   genContVar = D.genContinous++instance D.Entropy CauchyDistribution where+  entropy (CD _ s) = log s + log (4*pi)+    +instance D.MaybeEntropy CauchyDistribution where+  maybeEntropy = Just . D.entropy

@@ -21,7 +21,8 @@ import Data.Binary (Binary) import Data.Data (Data, Typeable) import GHC.Generics (Generic)-import Numeric.SpecFunctions (incompleteGamma,invIncompleteGamma,logGamma)+import Numeric.SpecFunctions (+  incompleteGamma,invIncompleteGamma,logGamma,digamma)  import qualified Statistics.Distribution         as D import qualified System.Random.MWC.Distributions as MWC@@ -68,6 +69,17 @@ instance D.MaybeVariance ChiSquared where     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance+    +instance D.Entropy ChiSquared where+  entropy (ChiSquared ndf) =+    let kHalf = 0.5 * fromIntegral ndf in+    kHalf +    + log 2 +    + logGamma kHalf+    + (1-kHalf) * digamma kHalf++instance D.MaybeEntropy ChiSquared where+  maybeEntropy = Just . D.entropy  instance D.ContGen ChiSquared where     genContVar (ChiSquared n) = MWC.chiSquare n

@@ -26,11 +26,13 @@ import Data.Binary (Binary) import Data.Data (Data, Typeable) import GHC.Generics (Generic)+import Numeric.MathFunctions.Constants (m_neg_inf) import qualified Statistics.Distribution         as D import qualified Statistics.Sample               as S import qualified System.Random.MWC.Distributions as MWC import Statistics.Types (Sample) + newtype ExponentialDistribution = ED {       edLambda :: Double     } deriving (Eq, Read, Show, Typeable, Data, Generic)@@ -42,7 +44,12 @@     complCumulative = complCumulative  instance D.ContDistr ExponentialDistribution where-    density  = density+    density (ED l) x+      | x < 0     = 0+      | otherwise = l * exp (-l * x)+    logDensity (ED l) x+      | x < 0     = m_neg_inf+      | otherwise = log l + (-l * x)     quantile = quantile  instance D.Mean ExponentialDistribution where@@ -60,6 +67,12 @@     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance +instance D.Entropy ExponentialDistribution where+  entropy (ED l) = 1 - log l++instance D.MaybeEntropy ExponentialDistribution where+  maybeEntropy = Just . D.entropy+ instance D.ContGen ExponentialDistribution where   genContVar = MWC.exponential . edLambda @@ -73,10 +86,6 @@                          | otherwise = exp (-l * x) {-# INLINE complCumulative #-} -density :: ExponentialDistribution -> Double -> Double-density (ED l) x | x < 0     = 0-                 | otherwise = l * exp (-l * x)-{-# INLINE density #-}  quantile :: ExponentialDistribution -> Double -> Double quantile (ED l) p

@@ -18,9 +18,11 @@  import Data.Binary (Binary) import Data.Data (Data, Typeable)+import Numeric.MathFunctions.Constants (m_neg_inf) import GHC.Generics (Generic) import qualified Statistics.Distribution as D-import Numeric.SpecFunctions (logBeta, incompleteBeta, invIncompleteBeta)+import Numeric.SpecFunctions (+  logBeta, incompleteBeta, invIncompleteBeta, digamma)   @@ -47,19 +49,23 @@   cumulative = cumulative  instance D.ContDistr FDistribution where-  density  = density+  density d x+    | x <= 0    = 0+    | otherwise = exp $ logDensity d x+  logDensity d x+    | x <= 0    = m_neg_inf+    | otherwise = logDensity d x   quantile = quantile  cumulative :: FDistribution -> Double -> Double cumulative (F n m _) x   | x <= 0       = 0   | isInfinite x = 1            -- Only matches +∞-  | x > 0        = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y))+  | otherwise    = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y)) -density :: FDistribution -> Double -> Double-density (F n m fac) x-  | x > 0     = exp $ fac + log x * (0.5 * n - 1) - log(m + n*x) * 0.5 * (n + m)-  | otherwise = 0+logDensity :: FDistribution -> Double -> Double+logDensity (F n m fac) x+  = fac + log x * (0.5 * n - 1) - log(m + n*x) * 0.5 * (n + m)  quantile :: FDistribution -> Double -> Double quantile (F n m _) p@@ -78,6 +84,19 @@   maybeStdDev (F n m _)     | m > 4     = Just $ 2 * sqr m * (m + n - 2) / (n * sqr (m - 2) * (m - 4))     | otherwise = Nothing++instance D.Entropy FDistribution where+  entropy (F n m _) =+    let nHalf = 0.5 * n+        mHalf = 0.5 * m in+    log (n/m) +    + logBeta nHalf mHalf+    + (1 - nHalf) * digamma nHalf +    - (1 + mHalf) * digamma mHalf+    + (nHalf + mHalf) * digamma (nHalf + mHalf)++instance D.MaybeEntropy FDistribution where+  maybeEntropy = Just . D.entropy  instance D.ContGen FDistribution where   genContVar = D.genContinous

@@ -28,8 +28,9 @@ import Data.Binary (Binary) import Data.Data (Data, Typeable) import GHC.Generics (Generic)-import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN)-import Numeric.SpecFunctions           (incompleteGamma, invIncompleteGamma)+import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf)+import Numeric.SpecFunctions (+  incompleteGamma, invIncompleteGamma, logGamma, digamma) import Statistics.Distribution.Poisson.Internal  as Poisson import qualified Statistics.Distribution         as D import qualified System.Random.MWC.Distributions as MWC@@ -67,6 +68,9 @@  instance D.ContDistr GammaDistribution where     density    = density+    logDensity (GD k theta) x+      | x <= 0    = m_neg_inf+      | otherwise = log x * (k - 1) - (x / theta) - logGamma k - log theta * k     quantile   = quantile  instance D.Variance GammaDistribution where@@ -83,6 +87,16 @@ instance D.MaybeVariance GammaDistribution where     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance++instance D.MaybeEntropy GammaDistribution where+  maybeEntropy (GD a l)+    | a > 0 && l > 0 = +      Just $+      a +      + log l +      + logGamma a +      + (1-a) * digamma a+    | otherwise = Nothing  instance D.ContGen GammaDistribution where     genContVar (GD a l) = MWC.gamma a l

@@ -8,29 +8,39 @@ -- Stability   : experimental -- Portability : portable ----- The Geometric distribution. This is the probability distribution of--- the number of Bernoulli trials needed to get one success, supported--- on the set [1,2..].+-- The Geometric distribution. There are two variants of+-- distribution. First is the probability distribution of the number+-- of Bernoulli trials needed to get one success, supported on the set+-- [1,2..] ('GeometricDistribution'). Sometimes it's referred to as+-- the /shifted/ geometric distribution to distinguish from another+-- one. ----- This distribution is sometimes referred to as the /shifted/--- geometric distribution, to distinguish it from a variant measuring--- the number of failures before the first success, defined over the--- set [0,1..].-+-- Second variant is probability distribution of the number of+-- failures before first success, defined over the set [0,1..]+-- ('GeometricDistribution0'). module Statistics.Distribution.Geometric     (       GeometricDistribution+    , GeometricDistribution0     -- * Constructors     , geometric+    , geometric0     -- ** Accessors     , gdSuccess+    , gdSuccess0     ) where -import Data.Binary (Binary)-import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Control.Monad  (liftM)+import Data.Binary    (Binary)+import Data.Data      (Data, Typeable)+import GHC.Generics   (Generic)+import Numeric.MathFunctions.Constants(m_pos_inf,m_neg_inf)+import qualified Statistics.Distribution         as D+import qualified System.Random.MWC.Distributions as MWC +----------------------------------------------------------------+-- Distribution over [1..]+ newtype GeometricDistribution = GD {       gdSuccess :: Double     } deriving (Eq, Read, Show, Typeable, Data, Generic)@@ -41,8 +51,14 @@     cumulative = cumulative  instance D.DiscreteDistr GeometricDistribution where-    probability = probability+    probability (GD s) n+      | n < 1     = 0+      | otherwise = s * (1-s) ** (fromIntegral n - 1)+    logProbability (GD s) n+       | n < 1     = m_neg_inf+       | otherwise = log s + log (1-s) * (fromIntegral n - 1) + instance D.Mean GeometricDistribution where     mean (GD s) = 1 / s     {-# INLINE mean #-}@@ -58,7 +74,22 @@     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance +instance D.Entropy GeometricDistribution where+  entropy (GD s)+    | s == 0 = m_pos_inf+    | s == 1 = 0+    | otherwise = negate $ (s * log s + (1-s) * log (1-s)) / s +instance D.MaybeEntropy GeometricDistribution where+  maybeEntropy = Just . D.entropy++instance D.DiscreteGen GeometricDistribution where+  genDiscreteVar (GD s) g = MWC.geometric1 s g+  {-# INLINE genDiscreteVar #-}+instance D.ContGen GeometricDistribution where+  genContVar d g = fromIntegral `liftM` D.genDiscreteVar d g+  {-# INLINE genContVar #-}+ -- | Create geometric distribution. geometric :: Double                -- ^ Success rate           -> GeometricDistribution@@ -68,11 +99,6 @@     error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x {-# INLINE geometric #-} -probability :: GeometricDistribution -> Int -> Double-probability (GD s) n | n < 1     = 0-                     | otherwise = s * (1-s) ** (fromIntegral n - 1)-{-# INLINE probability #-}- cumulative :: GeometricDistribution -> Double -> Double cumulative (GD s) x   | x < 1        = 0@@ -80,3 +106,57 @@   | isNaN      x = error "Statistics.Distribution.Geometric.cumulative: NaN input"   | otherwise    = 1 - (1-s) ^ (floor x :: Int) {-# INLINE cumulative #-}+++----------------------------------------------------------------+-- Distribution over [0..]++newtype GeometricDistribution0 = GD0 {+      gdSuccess0 :: Double+    } deriving (Eq, Read, Show, Typeable, Data, Generic)++instance Binary GeometricDistribution0++instance D.Distribution GeometricDistribution0 where+    cumulative (GD0 s) x = cumulative (GD s) (x + 1)++instance D.DiscreteDistr GeometricDistribution0 where+    probability    (GD0 s) n = D.probability    (GD s) (n + 1)+    logProbability (GD0 s) n = D.logProbability (GD s) (n + 1)++instance D.Mean GeometricDistribution0 where+    mean (GD0 s) = 1 / s - 1+    {-# INLINE mean #-}++instance D.Variance GeometricDistribution0 where+    variance (GD0 s) = D.variance (GD s)+    {-# INLINE variance #-}++instance D.MaybeMean GeometricDistribution0 where+    maybeMean = Just . D.mean++instance D.MaybeVariance GeometricDistribution0 where+    maybeStdDev   = Just . D.stdDev+    maybeVariance = Just . D.variance++instance D.Entropy GeometricDistribution0 where+  entropy (GD0 s) = D.entropy (GD s)++instance D.MaybeEntropy GeometricDistribution0 where+  maybeEntropy = Just . D.entropy++instance D.DiscreteGen GeometricDistribution0 where+  genDiscreteVar (GD0 s) g = MWC.geometric0 s g+  {-# INLINE genDiscreteVar #-}+instance D.ContGen GeometricDistribution0 where+  genContVar d g = fromIntegral `liftM` D.genDiscreteVar d g+  {-# INLINE genContVar #-}++-- | Create geometric distribution.+geometric0 :: Double                -- ^ Success rate+           -> GeometricDistribution0+geometric0 x+  | x >= 0 && x <= 1 = GD0 x+  | otherwise        =+    error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x+{-# INLINE geometric0 #-}

@@ -30,6 +30,7 @@ import Data.Binary (Binary) import Data.Data (Data, Typeable) import GHC.Generics (Generic)+import Numeric.MathFunctions.Constants (m_epsilon) import Numeric.SpecFunctions (choose) import qualified Statistics.Distribution as D @@ -60,7 +61,11 @@     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance -+instance D.Entropy HypergeometricDistribution where+  entropy = directEntropy+  +instance D.MaybeEntropy HypergeometricDistribution where+  maybeEntropy = Just . D.entropy  variance :: HypergeometricDistribution -> Double variance (HD m l k) = (k' * ml) * (1 - ml) * (l' - k') / (l' - 1)@@ -73,6 +78,14 @@ mean :: HypergeometricDistribution -> Double mean (HD m l k) = fromIntegral k * fromIntegral m / fromIntegral l {-# INLINE mean #-}++directEntropy :: HypergeometricDistribution -> Double+directEntropy d@(HD m _ _) =+    negate . sum $+  takeWhile (< negate m_epsilon) $+  dropWhile (not . (< negate m_epsilon)) $+  [ let x = probability d n in x * log x | n <- [0..m]]+  hypergeometric :: Int               -- ^ /m/                -> Int               -- ^ /l/

@@ -46,7 +46,7 @@     complCumulative = complCumulative  instance D.ContDistr NormalDistribution where-    density    = density+    logDensity = logDensity     quantile   = quantile  instance D.MaybeMean NormalDistribution where@@ -62,6 +62,12 @@ instance D.Variance NormalDistribution where     stdDev = stdDev +instance D.Entropy NormalDistribution where+  entropy d = 0.5 * log (2 * pi * exp 1 * D.variance d)++instance D.MaybeEntropy NormalDistribution where+  maybeEntropy = Just . D.entropy+ instance D.ContGen NormalDistribution where     genContVar d = MWC.normal (mean d) (stdDev d)     {-# INLINE genContVar #-}@@ -84,7 +90,7 @@ normalDistr m sd   | sd > 0    = ND { mean       = m                    , stdDev     = sd-                   , ndPdfDenom = m_sqrt_2_pi * sd+                   , ndPdfDenom = log $ m_sqrt_2_pi * sd                    , ndCdfDenom = m_sqrt_2 * sd                    }   | otherwise =@@ -99,8 +105,8 @@   where     (m,v) = S.meanVariance xs -density :: NormalDistribution -> Double -> Double-density d x = exp (-xm * xm / (2 * sd * sd)) / ndPdfDenom d+logDensity :: NormalDistribution -> Double -> Double+logDensity d x = (-xm * xm / (2 * sd * sd)) - ndPdfDenom d     where xm = x - mean d           sd = stdDev d 

@@ -29,8 +29,8 @@ import GHC.Generics (Generic) import qualified Statistics.Distribution as D import qualified Statistics.Distribution.Poisson.Internal as I-import Numeric.SpecFunctions (incompleteGamma)-+import Numeric.SpecFunctions (incompleteGamma,logFactorial)+import Numeric.MathFunctions.Constants (m_neg_inf)   newtype PoissonDistribution = PD {@@ -49,6 +49,9 @@  instance D.DiscreteDistr PoissonDistribution where     probability (PD lambda) x = I.probability lambda (fromIntegral x)+    logProbability (PD lambda) i+      | i < 0     = m_neg_inf+      | otherwise = log lambda * fromIntegral i - logFactorial i - lambda     {-# INLINE probability #-}  instance D.Variance PoissonDistribution where@@ -65,7 +68,12 @@ instance D.MaybeVariance PoissonDistribution where     maybeStdDev   = Just . D.stdDev +instance D.Entropy PoissonDistribution where+  entropy (PD lambda) = I.poissonEntropy lambda +instance D.MaybeEntropy PoissonDistribution where+  maybeEntropy = Just . D.entropy+ -- | Create Poisson distribution. poisson :: Double -> PoissonDistribution poisson l@@ -78,3 +86,6 @@ -- -- * Loader, C. (2000) Fast and Accurate Computation of Binomial --   Probabilities. <http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf>+-- * Adell, J., Lekuona, A., and Yu, Y. (2010) Sharp Bounds on the+--   Entropy of the Poisson Law and Related Quantities+--   <http://arxiv.org/pdf/1001.2897.pdf>

@@ -11,11 +11,12 @@  module Statistics.Distribution.Poisson.Internal     (-      probability+      probability, poissonEntropy     ) where -import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny)-import Numeric.SpecFunctions           (logGamma, stirlingError)+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@@ -31,3 +32,147 @@   | otherwise            = exp (-(stirlingError x) - bd0 x lambda) /                            (m_sqrt_2_pi * sqrt x) {-# INLINE probability #-}++-- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy+-- 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 - 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 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 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

@@ -21,7 +21,8 @@ import GHC.Generics (Generic) import qualified Statistics.Distribution as D import Statistics.Distribution.Transform (LinearTransform (..))-import Numeric.SpecFunctions (logBeta, incompleteBeta, invIncompleteBeta)+import Numeric.SpecFunctions (+  logBeta, incompleteBeta, invIncompleteBeta, digamma)  -- | Student-T distribution newtype StudentT = StudentT { studentTndf :: Double }@@ -39,7 +40,8 @@   cumulative = cumulative  instance D.ContDistr StudentT where-  density  = density+  density    d@(StudentT ndf) x = exp (logDensityUnscaled d x) / sqrt ndf+  logDensity d@(StudentT ndf) x = logDensityUnscaled d x - log (sqrt ndf)   quantile = quantile  cumulative :: StudentT -> Double -> Double@@ -49,9 +51,9 @@   where     ibeta = incompleteBeta (0.5 * ndf) 0.5 (ndf / (ndf + x*x)) -density :: StudentT -> Double -> Double-density (StudentT ndf) x =-    exp( log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf) ) / sqrt ndf+logDensityUnscaled :: StudentT -> Double -> Double+logDensityUnscaled (StudentT ndf) x =+    log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)  quantile :: StudentT -> Double -> Double quantile (StudentT ndf) p@@ -70,6 +72,15 @@ instance D.MaybeVariance StudentT where   maybeVariance (StudentT ndf) | ndf > 2   = Just $! ndf / (ndf - 2)                                | otherwise = Nothing++instance D.Entropy StudentT where+  entropy (StudentT ndf) =+    0.5 * (ndf+1) * (digamma ((1+ndf)/2) - digamma(ndf/2))+    + log (sqrt ndf) +    + logBeta (ndf/2) 0.5++instance D.MaybeEntropy StudentT where+  maybeEntropy = Just . D.entropy  instance D.ContGen StudentT where   genContVar = D.genContinous

@@ -55,7 +55,8 @@   cumulative (LinearTransform loc sc dist) x = D.cumulative dist $ (x-loc) / sc  instance D.ContDistr d => D.ContDistr (LinearTransform d) where-  density  (LinearTransform loc sc dist) x = D.density dist ((x-loc) / sc) / sc+  density    (LinearTransform loc sc dist) x = D.density    dist ((x-loc) / sc) / sc+  logDensity (LinearTransform loc sc dist) x = D.logDensity dist ((x-loc) / sc) - log sc   quantile (LinearTransform loc sc dist) p = loc + sc * D.quantile dist p  instance D.MaybeMean d => D.MaybeMean (LinearTransform d) where@@ -71,6 +72,14 @@ instance (D.Variance d) => D.Variance (LinearTransform d) where   variance (LinearTransform _ sc dist) = sc * sc * D.variance dist   stdDev   (LinearTransform _ sc dist) = sc * D.stdDev dist++instance (D.MaybeEntropy d, D.DiscreteDistr d) +         => D.MaybeEntropy (LinearTransform d) where+  maybeEntropy (LinearTransform _ _ dist) = D.maybeEntropy dist++instance (D.Entropy d, D.DiscreteDistr d) +         => D.Entropy (LinearTransform d) where+  entropy (LinearTransform _ _ dist) = D.entropy dist  instance D.ContGen d => D.ContGen (LinearTransform d) where   genContVar (LinearTransform loc sc d) g = do

@@ -73,5 +73,11 @@ instance D.MaybeVariance UniformDistribution where     maybeStdDev   = Just . D.stdDev +instance D.Entropy UniformDistribution where+  entropy (UniformDistribution a b) = log (b - a)++instance D.MaybeEntropy UniformDistribution where+  maybeEntropy = Just . D.entropy+ instance D.ContGen UniformDistribution where     genContVar (UniformDistribution a b) gen = MWC.uniformR (a,b) gen

@@ -203,7 +203,7 @@ data V = V {-# UNPACK #-} !Double {-# UNPACK #-} !Double  robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double-robustSumVar m samp = G.sum . G.map (square . subtract m) $ samp+robustSumVar m = G.sum . G.map (square . subtract m)   where square x = x * x {-# INLINE robustSumVar #-} 

@@ -1,5 +1,5 @@ name:           statistics-version:        0.10.4.1+version:        0.10.5.0 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful@@ -21,49 +21,6 @@   .   * Common statistical tests for significant differences between     samples.-  .-  Changes in 0.10.4.0-  .-  * Support for versions of GHC older than 7.2 is discontinued.-  .-  * All datatypes now support 'Data.Binary' and 'GHC.Generics'.-  .-  Changes in 0.10.3.0-  .-  * Bug fixes-  .-  Changes in 0.10.2.0-  .-  * Bugs in DCT and IDCT are fixed.-  .-  * Accesors for uniform distribution are added.-  .-  * 'ContGen' instances for all continous distribtuions are added.-  .-  * Beta distribution is added.-  .-  * Constructor for improper gamma distribtuion is added.-  .-  * Binomial distribution allows zero trials.-  .-  * Poisson distribution now accept zero parameter.-  .-  * Integer overflow in caculation of Wilcoxon-T test is fixed.-  .-  * Bug in 'ContGen' instance for normal distribution is fixed.-  .-  Changes in 0.10.1.0-  .-  * Kolmogorov-Smirnov nonparametric test added.-  .-  * Pearson's chi squared test added.-  .-  * Type class for generating random variates for given distribution-    is added.-  .-  * Modules 'Statistics.Math' and 'Statistics.Constants' are moved to-    the @math-functions@ package. They are still available but marked-    as deprecated.  license:        BSD3 license-file:   LICENSE@@ -137,7 +94,7 @@     deepseq >= 1.1.0.2,     erf,     monad-par         >= 0.3.4,-    mwc-random        >= 0.11.0.0,+    mwc-random        >= 0.13.0.0,     math-functions    >= 0.1.2,     primitive         >= 0.3,     vector            >= 0.7.1,

@@ -59,6 +59,7 @@    , discreteDistrTests (T :: T BinomialDistribution       )   , discreteDistrTests (T :: T GeometricDistribution      )+  , discreteDistrTests (T :: T GeometricDistribution0     )   , discreteDistrTests (T :: T HypergeometricDistribution )   , discreteDistrTests (T :: T PoissonDistribution        ) @@ -73,9 +74,10 @@ contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test contDistrTests t = testGroup ("Tests for: " ++ typeName t) $   cdfTests t ++-  [ testProperty "PDF sanity"              $ pdfSanityCheck   t-  , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t+  [ testProperty "PDF sanity"              $ pdfSanityCheck     t+  , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF   t   , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t+  , testProperty "log density check"       $ logDensityCheck    t   ]  -- Tests for discrete distribution@@ -85,6 +87,7 @@   [ testProperty "Prob. sanity"         $ probSanityCheck       t   , testProperty "CDF is sum of prob."  $ discreteCDFcorrect    t   , testProperty "Discrete CDF is OK"   $ cdfDiscreteIsCorrect  t+  , testProperty "log probabilty check" $ logProbabilityCheck   t   ]  -- Tests for distributions which have CDF@@ -144,7 +147,7 @@ -- CDF for discrete distribution uses <= for comparison cdfDiscreteIsCorrect :: (DiscreteDistr d) => T d -> d -> Property cdfDiscreteIsCorrect _ d-  = printTestCase (unlines $ map show badN)+  = printTestCase (unlines badN)   $ null badN     where     -- We are checking that:@@ -154,7 +157,7 @@     -- Apporixmate equality is tricky here. Scale is set by maximum     -- value of CDF and probability. Case when all proabilities are     -- zero should be trated specially.-    badN = [ (i,p,p1,dp, (p1-p-dp) / max p1 dp)+    badN = [ printf "N=%3i    p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp)            | i <- [0 .. 100]            , let p      = cumulative d $ fromIntegral i - 1e-6                  p1     = cumulative d $ fromIntegral i@@ -164,6 +167,19 @@            && relerr > 1e-14            ] +logDensityCheck :: (ContDistr d) => T d -> d -> Double -> Property+logDensityCheck _ d x+  = printTestCase (printf "density    = %g" p)+  $ printTestCase (printf "logDensity = %g" logP)+  $ printTestCase (printf "log p      = %g" (log p))+  $ printTestCase (printf "eps        = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))+  $ or [ p == 0     && logP == (-1/0)+       , p < 1e-308 && logP < 609+       , eq 1e-14 (log p) logP+       ]+  where+    p    = density d x+    logP = logDensity d x  -- PDF is positive pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool@@ -214,7 +230,21 @@     p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)     p2 = sum $ map (probability d) [n .. m] +logProbabilityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Property+logProbabilityCheck _ d x+  = printTestCase (printf "probability    = %g" p)+  $ printTestCase (printf "logProbability = %g" logP)+  $ printTestCase (printf "log p          = %g" (log p))+  $ printTestCase (printf "eps            = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))+  $ or [ p == 0     && logP == (-1/0)+       , p < 1e-308 && logP < 609+       , eq 1e-14 (log p) logP+       ]+  where+    p    = probability d x+    logP = logProbability d x +      ---------------------------------------------------------------- -- Arbitrary instances for ditributions@@ -230,6 +260,8 @@   arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10) instance QC.Arbitrary GeometricDistribution where   arbitrary = geometric <$> QC.choose (0,1)+instance QC.Arbitrary GeometricDistribution0 where+  arbitrary = geometric0 <$> QC.choose (0,1) instance QC.Arbitrary HypergeometricDistribution where   arbitrary = do l <- QC.choose (1,20)                  m <- QC.choose (0,l)