summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAlexeyKhudyakov <>2013-10-09 20:32:00 (GMT)
committerhdiff <hdiff@hdiff.luite.com>2013-10-09 20:32:00 (GMT)
commit11a18678d3a2e91c4ea0bfb2bc8a2a059622d52f (patch)
tree7d65ffcd31957ecf0047e89d39d6df2c5d7e4e6f
parentb71ca68286c816a79ace795ad745de5da06f4b96 (diff)
version 0.10.5.00.10.5.0
-rw-r--r--ChangeLog23
-rw-r--r--Statistics/Distribution.hs36
-rw-r--r--Statistics/Distribution/Beta.hs15
-rw-r--r--Statistics/Distribution/Binomial.hs25
-rw-r--r--Statistics/Distribution/CauchyLorentz.hs6
-rw-r--r--Statistics/Distribution/ChiSquared.hs14
-rw-r--r--Statistics/Distribution/Exponential.hs19
-rw-r--r--Statistics/Distribution/FDistribution.hs33
-rw-r--r--Statistics/Distribution/Gamma.hs18
-rw-r--r--Statistics/Distribution/Geometric.hs116
-rw-r--r--Statistics/Distribution/Hypergeometric.hs15
-rw-r--r--Statistics/Distribution/Normal.hs14
-rw-r--r--Statistics/Distribution/Poisson.hs15
-rw-r--r--Statistics/Distribution/Poisson/Internal.hs151
-rw-r--r--Statistics/Distribution/StudentT.hs21
-rw-r--r--Statistics/Distribution/Transform.hs11
-rw-r--r--Statistics/Distribution/Uniform.hs6
-rw-r--r--Statistics/Sample.hs2
-rw-r--r--statistics.cabal47
-rw-r--r--tests/Tests/Distribution.hs40
20 files changed, 520 insertions, 107 deletions
diff --git a/ChangeLog b/ChangeLog
index d27a474..dca54a5 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -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.
diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs
index 6809738..5c0894c 100644
--- a/Statistics/Distribution.hs
+++ b/Statistics/Distribution.hs
@@ -21,6 +21,8 @@ module Statistics.Distribution
, Mean(..)
, MaybeVariance(..)
, Variance(..)
+ , MaybeEntropy(..)
+ , Entropy(..)
-- ** Random number generation
, ContGen(..)
, DiscreteGen(..)
@@ -65,20 +67,36 @@ class Distribution d where
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+/&#948;/x/) equal to /density(x)/&#8901;&#948;/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/&#8804;/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 @@ class (Mean d, MaybeVariance d) => Variance d where
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.
diff --git a/Statistics/Distribution/Beta.hs b/Statistics/Distribution/Beta.hs
index 3c7245f..aff93b2 100644
--- a/Statistics/Distribution/Beta.hs
+++ b/Statistics/Distribution/Beta.hs
@@ -23,7 +23,8 @@ module Statistics.Distribution.Beta
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
@@ -82,6 +83,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
| a <= 0 || b <= 0 = m_NaN
diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs
index e8d522d..6827f4d 100644
--- a/Statistics/Distribution/Binomial.hs
+++ b/Statistics/Distribution/Binomial.hs
@@ -27,7 +27,9 @@ import Data.Binary (Binary)
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,6 +61,14 @@ instance D.MaybeVariance BinomialDistribution where
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
@@ -70,15 +80,13 @@ probability (BD n p) k
-- 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 #-}
@@ -90,6 +98,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
binomial :: Int -- ^ Number of trials.
diff --git a/Statistics/Distribution/CauchyLorentz.hs b/Statistics/Distribution/CauchyLorentz.hs
index a7b8129..bec1108 100644
--- a/Statistics/Distribution/CauchyLorentz.hs
+++ b/Statistics/Distribution/CauchyLorentz.hs
@@ -69,3 +69,9 @@ instance D.ContDistr CauchyDistribution where
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
diff --git a/Statistics/Distribution/ChiSquared.hs b/Statistics/Distribution/ChiSquared.hs
index c529747..1e87bd0 100644
--- a/Statistics/Distribution/ChiSquared.hs
+++ b/Statistics/Distribution/ChiSquared.hs
@@ -21,7 +21,8 @@ module Statistics.Distribution.ChiSquared (
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.MaybeMean ChiSquared where
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
diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs
index 01db6f2..3958d99 100644
--- a/Statistics/Distribution/Exponential.hs
+++ b/Statistics/Distribution/Exponential.hs
@@ -26,11 +26,13 @@ module Statistics.Distribution.Exponential
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 @@ instance D.Distribution ExponentialDistribution where
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 @@ instance D.MaybeVariance ExponentialDistribution where
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 @@ complCumulative (ED l) x | x <= 0 = 1
| 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
diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs
index 8dc2598..d3e3439 100644
--- a/Statistics/Distribution/FDistribution.hs
+++ b/Statistics/Distribution/FDistribution.hs
@@ -18,9 +18,11 @@ module Statistics.Distribution.FDistribution (
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 @@ instance D.Distribution FDistribution where
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
@@ -79,6 +85,19 @@ instance D.MaybeVariance FDistribution where
| 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
diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs
index b7af0ec..7e9507e 100644
--- a/Statistics/Distribution/Gamma.hs
+++ b/Statistics/Distribution/Gamma.hs
@@ -28,8 +28,9 @@ module Statistics.Distribution.Gamma
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.Distribution GammaDistribution where
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
@@ -84,6 +88,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
diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs
index 554b7bc..25e910d 100644
--- a/Statistics/Distribution/Geometric.hs
+++ b/Statistics/Distribution/Geometric.hs
@@ -8,28 +8,38 @@
-- 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
@@ -41,7 +51,13 @@ instance D.Distribution GeometricDistribution where
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
@@ -58,6 +74,21 @@ instance D.MaybeVariance GeometricDistribution where
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
@@ -68,11 +99,6 @@ geometric x
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 @@ cumulative (GD s) x
| 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 #-}
diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs
index b894b6b..a19466b 100644
--- a/Statistics/Distribution/Hypergeometric.hs
+++ b/Statistics/Distribution/Hypergeometric.hs
@@ -30,6 +30,7 @@ module Statistics.Distribution.Hypergeometric
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 @@ instance D.MaybeVariance HypergeometricDistribution where
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)
@@ -74,6 +79,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/
-> Int -- ^ /k/
diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs
index 1e3fb4b..10f763f 100644
--- a/Statistics/Distribution/Normal.hs
+++ b/Statistics/Distribution/Normal.hs
@@ -46,7 +46,7 @@ instance D.Distribution NormalDistribution where
complCumulative = complCumulative
instance D.ContDistr NormalDistribution where
- density = density
+ logDensity = logDensity
quantile = quantile
instance D.MaybeMean NormalDistribution where
@@ -62,6 +62,12 @@ instance D.MaybeVariance NormalDistribution where
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 :: Double -- ^ Mean of distribution
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 @@ normalFromSample xs
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
diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs
index 96f713b..87f822e 100644
--- a/Statistics/Distribution/Poisson.hs
+++ b/Statistics/Distribution/Poisson.hs
@@ -29,8 +29,8 @@ import Data.Data (Data, Typeable)
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.Distribution PoissonDistribution where
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,6 +68,11 @@ instance D.MaybeMean PoissonDistribution where
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
@@ -78,3 +86,6 @@ poisson l
--
-- * 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>
diff --git a/Statistics/Distribution/Poisson/Internal.hs b/Statistics/Distribution/Poisson/Internal.hs
index 91f5a4b..039adbf 100644
--- a/Statistics/Distribution/Poisson/Internal.hs
+++ b/Statistics/Distribution/Poisson/Internal.hs
@@ -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 @@ probability lambda x
| 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 \ No newline at end of file
diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs
index 1439371..03fcaba 100644
--- a/Statistics/Distribution/StudentT.hs
+++ b/Statistics/Distribution/StudentT.hs
@@ -21,7 +21,8 @@ import Data.Data (Data, Typeable)
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 @@ instance D.Distribution StudentT where
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 @@ cumulative (StudentT ndf) x
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
@@ -71,6 +73,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
diff --git a/Statistics/Distribution/Transform.hs b/Statistics/Distribution/Transform.hs
index ee72e4e..f575091 100644
--- a/Statistics/Distribution/Transform.hs
+++ b/Statistics/Distribution/Transform.hs
@@ -55,7 +55,8 @@ instance D.Distribution d => D.Distribution (LinearTransform d) where
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
@@ -72,6 +73,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
x <- D.genContVar d g
diff --git a/Statistics/Distribution/Uniform.hs b/Statistics/Distribution/Uniform.hs
index eb7d018..95a2b79 100644
--- a/Statistics/Distribution/Uniform.hs
+++ b/Statistics/Distribution/Uniform.hs
@@ -73,5 +73,11 @@ instance D.MaybeMean UniformDistribution where
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
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
index ab98006..a068f60 100644
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -203,7 +203,7 @@ kurtosis xs = c4 / (c2 * c2) - 3
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 #-}
diff --git a/statistics.cabal b/statistics.cabal
index 6f22d17..665f3f6 100644
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -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 @@ description:
.
* 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 @@ library
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,
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
index ad2a1f3..bfc2ba7 100644
--- a/tests/Tests/Distribution.hs
+++ b/tests/Tests/Distribution.hs
@@ -59,6 +59,7 @@ distributionTests = testGroup "Tests for all distributions"
, discreteDistrTests (T :: T BinomialDistribution )
, discreteDistrTests (T :: T GeometricDistribution )
+ , discreteDistrTests (T :: T GeometricDistribution0 )
, discreteDistrTests (T :: T HypergeometricDistribution )
, discreteDistrTests (T :: T PoissonDistribution )
@@ -73,9 +74,10 @@ distributionTests = testGroup "Tests for all distributions"
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 @@ discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
[ 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 @@ cdfComplementIsCorrect _ d x = (eq 1e-14) 1 (cumulative d x + complCumulative d
-- 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 @@ cdfDiscreteIsCorrect _ d
-- 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 @@ cdfDiscreteIsCorrect _ d
&& 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,6 +230,20 @@ discreteCDFcorrect _ d a b
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
+
----------------------------------------------------------------
@@ -230,6 +260,8 @@ instance QC.Arbitrary BetaDistribution where
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)