diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs --- a/Statistics/Distribution.hs +++ b/Statistics/Distribution.hs @@ -8,36 +8,60 @@ -- Stability : experimental -- Portability : portable -- --- Types and functions common to many probability distributions. +-- Types classes for probability distrubutions module Statistics.Distribution ( + -- * Type classes Distribution(..) + , DiscreteDistr(..) + , ContDistr(..) , Mean(..) , Variance(..) + -- * Helper functions , findRoot + , sumProbabilities ) where --- | The interface shared by all probability distributions. -class Distribution d where - -- | Probability density function. The probability that a - -- the random variable /X/ has the value /x/, i.e. P(/X/=/x/). - density :: d -> Double -> Double +import qualified Data.Vector.Unboxed as U +-- | Type class common to all distributions. Only c.d.f. could be +-- defined for both discrete and continous distributions. +class Distribution d where -- | Cumulative distribution function. The probability that a - -- random variable /X/ is less than /x/, i.e. P(/X/≤/x/). + -- random variable /X/ is less or equal than /x/, + -- i.e. P(/X/≤/x/). cumulative :: d -> Double -> Double - -- | Inverse of the cumulative distribution function. The value - -- /x/ for which P(/X/≤/x/). + +-- | Discrete probability distribution. +class Distribution d => DiscreteDistr d where + -- | Probability of n-th outcome. + probability :: d -> Int -> Double + + +-- | Continuous probability distributuion +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 + + -- | Inverse of the cumulative distribution function. The value + -- /x/ for which P(/X/≤/x/) = /p/. quantile :: d -> Double -> Double + +-- | Type class for distributions with mean. class Distribution d => Mean d where mean :: d -> Double + +-- | Type class for distributions with variance. class Mean d => Variance d where variance :: d -> Double + data P = P {-# UNPACK #-} !Double {-# UNPACK #-} !Double -- | Approximate the value of /X/ for which P(/x/>/X/)=/p/. @@ -46,7 +70,8 @@ -- bisection with the given guess as a starting point. The upper and -- lower bounds specify the interval in which the probability -- distribution reaches the value /p/. -findRoot :: Distribution d => d +findRoot :: ContDistr d => + d -- ^ Distribution -> Double -- ^ Probability /p/ -> Double -- ^ Initial guess -> Double -- ^ Lower bound on interval @@ -70,3 +95,11 @@ | otherwise = P dx' x' accuracy = 1e-15 maxIters = 150 + +-- | Sum probabilities in inclusive interval. +sumProbabilities :: DiscreteDistr d => d -> Int -> Int -> Double +sumProbabilities d low hi = + -- Return value is forced to be less than 1 to guard againist roundoff errors. + -- ATTENTION! this check should be removed for testing or it could mask bugs. + min 1 . U.sum . U.map (probability d) $ U.enumFromTo low hi +{-# INLINE sumProbabilities #-} diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs --- a/Statistics/Distribution/Binomial.hs +++ b/Statistics/Distribution/Binomial.hs @@ -24,13 +24,9 @@ ) where import Control.Exception (assert) -import qualified Data.Vector.Unboxed as U -import Data.Int (Int64) import Data.Typeable (Typeable) -import Statistics.Constants (m_epsilon) import qualified Statistics.Distribution as D -import Statistics.Distribution.Normal (standard) -import Statistics.Math (choose, logFactorial) +import Statistics.Math (choose) -- | The binomial distribution. data BinomialDistribution = BD { @@ -41,80 +37,37 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution BinomialDistribution where - density = density cumulative = cumulative - quantile = quantile +instance D.DiscreteDistr BinomialDistribution where + probability = probability + instance D.Variance BinomialDistribution where variance = variance instance D.Mean BinomialDistribution where mean = mean -density :: BinomialDistribution -> Double -> Double -density (BD n p) x - | not (isIntegral x) = integralError "density" - | n == 0 = 1 - | x < 0 || x > n' = 0 - | n <= 50 || x < 2 = sign * p'' ** x' * (n `choose` fx) * q'' ** nx' - | otherwise = sign * exp (x' * log p'' + nx' * log q'' + lf) - where sign = oddX * oddNX - (x',p',q') | x > n' / 2 = (n'-x, q, p) - | otherwise = (x, p, q) - oddX | p' < 0 && odd fx = -1 - | otherwise = 1 - oddNX | q' < 0 && odd nx = -1 - | otherwise = 1 - p'' = abs p' - q'' = abs q' - q = 1 - p - nx = n - fx - nx' = fromIntegral nx - fx = floor x' - n' = fromIntegral n - lf = logFactorial n - logFactorial nx - logFactorial fx -cumulative :: BinomialDistribution -> Double -> Double -cumulative d x - | isIntegral x = U.sum . U.map (density d . fromIntegral) . U.enumFromTo (0::Int) . floor $ x - | otherwise = integralError "cumulative" - -isIntegral :: Double -> Bool -isIntegral x = x == floorf x - -floorf :: Double -> Double -floorf = fromIntegral . (floor :: Double -> Int64) +-- This could be slow for bin n +probability :: BinomialDistribution -> Int -> Double +probability (BD n p) k + | k < 0 || k > n = 0 + | n == 0 = 1 + | otherwise = choose n k * p^k * (1-p)^(n-k) +{-# INLINE probability #-} -quantile :: BinomialDistribution -> Double -> Double -quantile dist@(BD n p) prob - | isNaN prob = prob - | p == 1 = n' - | n' < 1e5 = fst (search 1 y0 z0) - | otherwise = let dy = floorf (n' / 1000) - in narrow dy (search dy y0 z0) - where q = 1 - p - n' = fromIntegral n - y0 = n' `min` floorf (µ + σ * (d + γ * (d * d - 1) / 6) + 0.5) - where µ = n' * p - σ = sqrt (n' * p * q) - d = D.quantile standard prob - γ = (q - p) / σ - z0 = cumulative dist y0 - search dy y1 z1 | z0 >= prob' = left y1 z1 - | otherwise = right y1 - where - prob' = prob * (1 - 64 * m_epsilon) - left y oldZ | y == 0 || z < prob' = (y, oldZ) - | otherwise = left (max 0 y') z - where z = cumulative dist y' - y' = y - dy - right y | y' >= n' || z >= prob' = (y', z) - | otherwise = right y' - where z = cumulative dist y' - y' = y + dy - narrow dy (y,z) | dy <= 1 || dy' <= n'/1e15 = y - | otherwise = narrow dy' (search dy y z) - where dy' = floorf (dy / 100) +-- Summation from different sides required to reduce roundoff errors +cumulative :: BinomialDistribution -> Double -> Double +cumulative d@(BD n _) x + | k < 0 = 0 + | k >= n = 1 + | k < m = D.sumProbabilities d 0 k + | otherwise = 1 - D.sumProbabilities d (k+1) n + where + m = floor (mean d) + k = floor x +{-# INLINE cumulative #-} mean :: BinomialDistribution -> Double mean (BD n p) = fromIntegral n * p @@ -124,6 +77,7 @@ variance (BD n p) = fromIntegral n * p * (1 - p) {-# INLINE variance #-} +-- | Construct binomial distribution binomial :: Int -- ^ Number of trials. -> Double -- ^ Probability. -> BinomialDistribution @@ -132,7 +86,3 @@ assert (p > 0 && p < 1) $ BD n p {-# INLINE binomial #-} - -integralError :: String -> a -integralError f = error ("Statistics.Distribution.Binomial." ++ f ++ - ": non-integer-valued input") diff --git a/Statistics/Distribution/ChiSquared.hs b/Statistics/Distribution/ChiSquared.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/ChiSquared.hs @@ -0,0 +1,80 @@ +{-# LANGUAGE DeriveDataTypeable #-} +-- | +-- Module : Statistics.Distribution.ChiSquared +-- Copyright : (c) 2010 Alexey Khudyakov +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- The chi-squared distribution. This is a continuous probability +-- distribution of sum of squares of k independent standard normal +-- distributions. It's commonly used in statistical tests +module Statistics.Distribution.ChiSquared ( + ChiSquared + -- Constructors + , chiSquared + , chiSquaredNDF + ) where + +import Data.Typeable (Typeable) +import Statistics.Constants (m_huge) +import Statistics.Math (incompleteGamma,logGamma) + +import qualified Statistics.Distribution as D + + +-- | Chi-squared distribution +newtype ChiSquared = ChiSquared Int + deriving (Show,Typeable) + +-- | Get number of degrees of freedom +chiSquaredNDF :: ChiSquared -> Int +chiSquaredNDF (ChiSquared ndf) = ndf +{-# INLINE chiSquaredNDF #-} + +-- | Construct chi-squared distribution. Number of degrees of free +chiSquared :: Int -> ChiSquared +chiSquared x = ChiSquared x +{-# INLINE chiSquared #-} + +instance D.Distribution ChiSquared where + cumulative = cumulative + +instance D.ContDistr ChiSquared where + density = density + quantile = quantile + +instance D.Mean ChiSquared where + mean (ChiSquared ndf) = fromIntegral ndf + {-# INLINE mean #-} + +instance D.Variance ChiSquared where + variance (ChiSquared ndf) = fromIntegral (2*ndf) + {-# INLINE variance #-} + +cumulative :: ChiSquared -> Double -> Double +cumulative chi x + | x <= 0 = 0 + | otherwise = incompleteGamma (ndf/2) (x/2) + where + ndf = fromIntegral $ chiSquaredNDF chi +{-# INLINE cumulative #-} + +density :: ChiSquared -> Double -> Double +density chi x + | x <= 0 = 0 + | otherwise = exp $ log x * (ndf2 - 1) - x2 - logGamma ndf2 - log 2 * ndf2 + where + ndf = fromIntegral $ chiSquaredNDF chi + ndf2 = ndf/2 + x2 = x/2 +{-# INLINE density #-} + +quantile :: ChiSquared -> Double -> Double +quantile d@(ChiSquared ndf) p + | p == 0 = -1/0 + | p == 1 = 1/0 + | otherwise = D.findRoot d p (fromIntegral ndf) 0 m_huge +{-# INLINE quantile #-} diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs --- a/Statistics/Distribution/Exponential.hs +++ b/Statistics/Distribution/Exponential.hs @@ -17,8 +17,8 @@ ( ExponentialDistribution -- * Constructors - , fromLambda - , fromSample + , exponential + , exponentialFromSample -- * Accessors , edLambda ) where @@ -33,13 +33,12 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution ExponentialDistribution where - density (ED l) x = l * exp (-l * x) - {-# INLINE density #-} - cumulative (ED l) x = 1 - exp (-l * x) - {-# INLINE cumulative #-} - quantile (ED l) p = -log (1 - p) / l - {-# INLINE quantile #-} + cumulative = cumulative +instance D.ContDistr ExponentialDistribution where + density = density + quantile = quantile + instance D.Variance ExponentialDistribution where variance (ED l) = 1 / (l * l) {-# INLINE variance #-} @@ -48,11 +47,28 @@ mean (ED l) = 1 / l {-# INLINE mean #-} -fromLambda :: Double -- ^ λ (scale) parameter. - -> ExponentialDistribution -fromLambda = ED -{-# INLINE fromLambda #-} +cumulative :: ExponentialDistribution -> Double -> Double +cumulative (ED l) x | x < 0 = 0 + | otherwise = 1 - exp (-l * x) +{-# INLINE cumulative #-} -fromSample :: Sample -> ExponentialDistribution -fromSample = ED . S.mean -{-# INLINE fromSample #-} +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 = -log (1 - p) / l +{-# INLINE quantile #-} + +-- | Create exponential distribution +exponential :: Double -- ^ λ (scale) parameter. + -> ExponentialDistribution +exponential = ED +{-# INLINE exponential #-} + +-- | Create exponential distribution from sample. No tests are made to +-- check whether it really exponential +exponentialFromSample :: Sample -> ExponentialDistribution +exponentialFromSample = ED . S.mean +{-# INLINE exponentialFromSample #-} diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs --- a/Statistics/Distribution/Gamma.hs +++ b/Statistics/Distribution/Gamma.hs @@ -18,9 +18,7 @@ ( GammaDistribution -- * Constructors - --, fromParams - --, fromSample - --, standard + , gammaDistr -- * Accessors , gdShape , gdScale @@ -37,9 +35,15 @@ , gdScale :: {-# UNPACK #-} !Double -- ^ Scale parameter, ϑ. } deriving (Eq, Read, Show, Typeable) +gammaDistr :: Double -> Double -> GammaDistribution +gammaDistr = GD +{-# INLINE gammaDistr #-} + instance D.Distribution GammaDistribution where - density = density cumulative = cumulative + +instance D.ContDistr GammaDistribution where + density = density quantile = quantile instance D.Variance GammaDistribution where @@ -51,11 +55,15 @@ {-# INLINE mean #-} density :: GammaDistribution -> Double -> Double -density (GD a l) x = x ** (a-1) * exp (-x/l) / (exp (logGamma a) * l ** a) +density (GD a l) x + | x <= 0 = 0 + | otherwise = x ** (a-1) * exp (-x/l) / (exp (logGamma a) * l ** a) {-# INLINE density #-} cumulative :: GammaDistribution -> Double -> Double -cumulative (GD a l) x = incompleteGamma a (x/l) / exp (logGamma a) +cumulative (GD k l) x + | x <= 0 = 0 + | otherwise = incompleteGamma k (x/l) {-# INLINE cumulative #-} quantile :: GammaDistribution -> Double -> Double diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs --- a/Statistics/Distribution/Geometric.hs +++ b/Statistics/Distribution/Geometric.hs @@ -21,9 +21,9 @@ ( GeometricDistribution -- * Constructors - , fromSuccess + , geometric -- ** Accessors - , pdSuccess + , gdSuccess ) where import Control.Exception (assert) @@ -31,14 +31,15 @@ import qualified Statistics.Distribution as D newtype GeometricDistribution = GD { - pdSuccess :: Double + gdSuccess :: Double } deriving (Eq, Read, Show, Typeable) instance D.Distribution GeometricDistribution where - density = density cumulative = cumulative - quantile = quantile +instance D.DiscreteDistr GeometricDistribution where + probability = probability + instance D.Variance GeometricDistribution where variance (GD s) = (1 - s) / (s * s) {-# INLINE variance #-} @@ -47,19 +48,19 @@ mean (GD s) = 1 / s {-# INLINE mean #-} -fromSuccess :: Double -> GeometricDistribution -fromSuccess x = assert (x >= 0 && x <= 1) - GD x -{-# INLINE fromSuccess #-} +-- | Create geometric distribution +geometric :: Double -- ^ Success rate + -> GeometricDistribution +geometric x = assert (x >= 0 && x <= 1) + GD x +{-# INLINE geometric #-} -density :: GeometricDistribution -> Double -> Double -density (GD s) x = s * (1-s) ** (x-1) -{-# INLINE density #-} +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 = 1 - (1-s) ** x +cumulative (GD s) x | x < 1 = 0 + | otherwise = 1 - (1-s) ^ (floor x :: Int) {-# INLINE cumulative #-} - -quantile :: GeometricDistribution -> Double -> Double -quantile (GD s) p = log (1 - p) / log (1 - s) -{-# INLINE quantile #-} diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs --- a/Statistics/Distribution/Hypergeometric.hs +++ b/Statistics/Distribution/Hypergeometric.hs @@ -20,7 +20,7 @@ ( HypergeometricDistribution -- * Constructors - , fromParams + , hypergeometric -- ** Accessors , hdM , hdL @@ -28,10 +28,8 @@ ) where import Control.Exception (assert) -import qualified Data.Vector.Unboxed as U import Data.Typeable (Typeable) -import Statistics.Math (choose, logFactorial) -import Statistics.Constants (m_max_exp) +import Statistics.Math (choose) import qualified Statistics.Distribution as D data HypergeometricDistribution = HD { @@ -41,10 +39,11 @@ } deriving (Eq, Read, Show, Typeable) instance D.Distribution HypergeometricDistribution where - density = density - cumulative = cumulative - quantile = quantile + cumulative d x = D.sumProbabilities d 0 (floor x) +instance D.DiscreteDistr HypergeometricDistribution where + probability = probability + instance D.Variance HypergeometricDistribution where variance = variance @@ -63,45 +62,21 @@ mean (HD m l k) = fromIntegral k * fromIntegral m / fromIntegral l {-# INLINE mean #-} -fromParams :: Int -- ^ /m/ - -> Int -- ^ /l/ - -> Int -- ^ /k/ - -> HypergeometricDistribution -fromParams m l k = - assert (m > 0 && m <= l) . +hypergeometric :: Int -- ^ /m/ + -> Int -- ^ /l/ + -> Int -- ^ /k/ + -> HypergeometricDistribution +hypergeometric m l k = + assert (m >= 0 && m <= l) . assert (l > 0) . assert (k > 0 && k <= l) $ HD m l k -{-# INLINE fromParams #-} - -density :: HypergeometricDistribution -> Double -> Double -density (HD mi li ki) x - | l <= 70 = (mi <> xi) * ((li - mi) <> (ki - xi)) / (li <> ki) - | r > maxVal = 1/0 - | otherwise = exp r - where - a <> b = a `choose` b - r = f m + f (l-m) - f l - f xi - f (k-xi) + f k - - f (m-xi) - f (l-m-k+xi) + f (l-k) - f = logFactorial - maxVal = fromIntegral (m_max_exp - 1) * log 2 - xi = floor x - m = fromIntegral mi - l = fromIntegral li - k = fromIntegral ki -{-# INLINE density #-} - -cumulative :: HypergeometricDistribution -> Double -> Double -cumulative d@(HD m l k) x - | x < fromIntegral imin = 0 - | x >= fromIntegral imax = 1 - | otherwise = min r 1 - where - imin = max 0 (k - l + m) - imax = min k m - r = U.sum . U.map (density d . fromIntegral) . U.enumFromTo imin . floor $ x -{-# INLINE cumulative #-} +{-# INLINE hypergeometric #-} -quantile :: HypergeometricDistribution -> Double -> Double -quantile = error "Statistics.Distribution.Hypergeometric.quantile: not yet implemented" -{-# INLINE quantile #-} +-- Naive implementation +probability :: HypergeometricDistribution -> Int -> Double +probability (HD mi li ki) n + | n < max 0 (mi+ki-li) || n > min mi ki = 0 + | otherwise = + choose mi n * choose (li - mi) (ki - n) / choose li ki +{-# INLINE probability #-} diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs --- a/Statistics/Distribution/Normal.hs +++ b/Statistics/Distribution/Normal.hs @@ -15,8 +15,8 @@ ( NormalDistribution -- * Constructors - , fromParams - , fromSample + , normalDistr + , normalFromSample , standard ) where @@ -29,15 +29,17 @@ -- | The normal distribution. data NormalDistribution = ND { - mean :: {-# UNPACK #-} !Double - , variance :: {-# UNPACK #-} !Double + mean :: {-# UNPACK #-} !Double + , variance :: {-# UNPACK #-} !Double , ndPdfDenom :: {-# UNPACK #-} !Double , ndCdfDenom :: {-# UNPACK #-} !Double } deriving (Eq, Read, Show, Typeable) instance D.Distribution NormalDistribution where - density = density cumulative = cumulative + +instance D.ContDistr NormalDistribution where + density = density quantile = quantile instance D.Variance NormalDistribution where @@ -48,38 +50,36 @@ -- | Standard normal distribution with mean equal to 0 and variance equal to 1 standard :: NormalDistribution -standard = ND { - mean = 0.0 - , variance = 1.0 - , ndPdfDenom = m_sqrt_2_pi - , ndCdfDenom = m_sqrt_2 - } +standard = ND { mean = 0.0 + , variance = 1.0 + , ndPdfDenom = m_sqrt_2_pi + , ndCdfDenom = m_sqrt_2 + } -- | Create normal distribution from parameters -fromParams :: Double -- ^ Mean of distribution - -> Double -- ^ Variance of distribution - -> NormalDistribution -fromParams m v = assert (v > 0) - ND { - mean = m - , variance = v - , ndPdfDenom = m_sqrt_2_pi * sv - , ndCdfDenom = m_sqrt_2 * sv - } +normalDistr :: Double -- ^ Mean of distribution + -> Double -- ^ Variance of distribution + -> NormalDistribution +normalDistr m v = assert (v > 0) + ND { mean = m + , variance = v + , ndPdfDenom = m_sqrt_2_pi * sv + , ndCdfDenom = m_sqrt_2 * sv + } where sv = sqrt v -- | Create distribution using parameters estimated from -- sample. Variance is estimated using maximum likelihood method -- (biased estimation). -fromSample :: S.Sample -> NormalDistribution -fromSample a = fromParams (S.mean a) (S.variance a) +normalFromSample :: S.Sample -> NormalDistribution +normalFromSample a = normalDistr (S.mean a) (S.variance a) density :: NormalDistribution -> Double -> Double density d x = exp (-xm * xm / (2 * variance d)) / ndPdfDenom d where xm = x - mean d cumulative :: NormalDistribution -> Double -> Double -cumulative d x = erfc (-(x-mean d) / ndCdfDenom d) / 2 +cumulative d x = erfc ((mean d - x) / ndCdfDenom d) / 2 quantile :: NormalDistribution -> Double -> Double quantile d p diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs --- a/Statistics/Distribution/Poisson.hs +++ b/Statistics/Distribution/Poisson.hs @@ -17,53 +17,44 @@ ( PoissonDistribution -- * Constructors - , fromLambda - -- , fromSample + , poisson + -- * Accessors + , poissonLambda ) where import Data.Typeable (Typeable) -import qualified Data.Vector.Unboxed as U import qualified Statistics.Distribution as D -import Statistics.Constants (m_huge) -import Statistics.Math (factorial, logGamma) +import Statistics.Math (logGamma, factorial) newtype PoissonDistribution = PD { - pdLambda :: Double + poissonLambda :: Double } deriving (Eq, Read, Show, Typeable) instance D.Distribution PoissonDistribution where - density = density - cumulative = cumulative - quantile = quantile + cumulative d x = D.sumProbabilities d 0 (floor x) + {-# INLINE cumulative #-} +instance D.DiscreteDistr PoissonDistribution where + probability = probability + instance D.Variance PoissonDistribution where - variance = pdLambda + variance = poissonLambda {-# INLINE variance #-} instance D.Mean PoissonDistribution where - mean = pdLambda + mean = poissonLambda {-# INLINE mean #-} -fromLambda :: Double -> PoissonDistribution -fromLambda = PD -{-# INLINE fromLambda #-} - -density :: PoissonDistribution -> Double -> Double -density (PD l) x - | x < 0 = 0 - | l >= 100 && x >= l * 10 = 0 - | l >= 3 && x >= l * 100 = 0 - | x >= max 1 l * 200 = 0 - | l < 20 && x <= 100 = exp (-l) * l ** x / factorial (floor x) - | otherwise = exp (x * log l - logGamma (x + 1) - l) -{-# INLINE density #-} - -cumulative :: PoissonDistribution -> Double -> Double -cumulative d = U.sum . U.map (density d . fromIntegral) . - U.enumFromTo (0::Int) . floor -{-# INLINE cumulative #-} +-- | Create po +poisson :: Double -> PoissonDistribution +poisson = PD +{-# INLINE poisson #-} -quantile :: PoissonDistribution -> Double -> Double -quantile d p = fromIntegral . r $ D.findRoot d p (pdLambda d) 0 m_huge - where r = round :: Double -> Int -{-# INLINE quantile #-} +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 #-} diff --git a/Statistics/Math.hs b/Statistics/Math.hs --- a/Statistics/Math.hs +++ b/Statistics/Math.hs @@ -88,29 +88,26 @@ -- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For -- values of /k/ > 30, this uses an approximation for performance --- reasons. The approximation is accurate to 7 decimal places in the --- worst case, but is typically accurate to 9 decimal places or --- better. +-- reasons. The approximation is accurate to 12 decimal places in the +-- worst case -- -- Example: -- -- > 7 `choose` 3 == 35 choose :: Int -> Int -> Double n `choose` k - | k > n = 0 - | k < 30 = U.foldl' go 1 . U.enumFromTo 1 $ k' + | k > n = 0 + | k' < 50 = U.foldl' go 1 . U.enumFromTo 1 $ k' | approx < max64 = fromIntegral . round64 $ approx | otherwise = approx where - approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k) + k' = min k (n-k) + approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k') -- Less numerically stable: -- exp $ lg (n+1) - lg (k+1) - lg (n-k+1) -- where lg = logGamma . fromIntegral go a i = a * (nk + j) / j where j = fromIntegral i :: Double - k' | n_k < k = n_k - | otherwise = k - where n_k = n - k nk = fromIntegral (n - k') max64 = fromIntegral (maxBound :: Int64) round64 x = round x :: Int64 @@ -145,11 +142,11 @@ -- | Compute the normalized lower incomplete gamma function -- γ(/s/,/x/). Normalization means that --- γ(∞,/x/)=1. Uses Algorithm AS 239 by Shea. +-- γ(/s/,∞)=1. Uses Algorithm AS 239 by Shea. incompleteGamma :: Double -- ^ /s/ -> Double -- ^ /x/ -> Double -incompleteGamma x p +incompleteGamma p x | x < 0 || p <= 0 = m_pos_inf | x == 0 = 0 | p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1)) diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs --- a/Statistics/Sample.hs +++ b/Statistics/Sample.hs @@ -38,6 +38,8 @@ -- $robust , variance , varianceUnbiased + , meanVariance + , meanVarianceUnb , stdDev , varianceWeighted @@ -204,17 +206,15 @@ sqr :: Double -> Double sqr x = x * x -robustSumVar :: (G.Vector v Double) => v Double -> Double -robustSumVar samp = G.sum . G.map (sqr . subtract m) $ samp - where - m = mean samp +robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double +robustSumVar m samp = G.sum . G.map (sqr . subtract m) $ samp {-# INLINE robustSumVar #-} -- | Maximum likelihood estimate of a sample's variance. Also known -- as the population variance, where the denominator is /n/. variance :: (G.Vector v Double) => v Double -> Double variance samp - | n > 1 = robustSumVar samp / fromIntegral n + | n > 1 = robustSumVar (mean samp) samp / fromIntegral n | otherwise = 0 where n = G.length samp @@ -224,11 +224,35 @@ -- sample variance, where the denominator is /n/-1. varianceUnbiased :: (G.Vector v Double) => v Double -> Double varianceUnbiased samp - | n > 1 = robustSumVar samp / fromIntegral (n-1) + | n > 1 = robustSumVar (mean samp) samp / fromIntegral (n-1) | otherwise = 0 where n = G.length samp {-# INLINE varianceUnbiased #-} + +-- | Calculate mean and maximum likelihood estimate of variance. This +-- function should be used if both mean and variance are required +-- since it will calculate mean only once. +meanVariance :: (G.Vector v Double) => v Double -> (Double,Double) +meanVariance samp + | n > 1 = (m, robustSumVar m samp / fromIntegral n) + | otherwise = (m, 0) + where + n = G.length samp + m = mean samp +{-# INLINE meanVariance #-} + +-- | Calculate mean and unbiased estimate of variance. This +-- function should be used if both mean and variance are required +-- since it will calculate mean only once. +meanVarianceUnb :: (G.Vector v Double) => v Double -> (Double,Double) +meanVarianceUnb samp + | n > 1 = (m, robustSumVar m samp / fromIntegral (n-1)) + | otherwise = (m, 0) + where + n = G.length samp + m = mean samp +{-# INLINE meanVarianceUnb #-} -- | Standard deviation. This is simply the square root of the -- unbiased estimate of the variance. diff --git a/Statistics/Test/NonParametric.hs b/Statistics/Test/NonParametric.hs new file mode 100644 --- /dev/null +++ b/Statistics/Test/NonParametric.hs @@ -0,0 +1,325 @@ +-- | +-- Module : Statistics.Test.NonParametric +-- Copyright : (c) 2010 Neil Brown +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Functions for performing non-parametric tests (i.e. tests without an assumption +-- of underlying distribution). +module Statistics.Test.NonParametric + (-- * Mann-Whitney U test (non-parametric equivalent to the independent t-test) + mannWhitneyU, mannWhitneyUCriticalValue, mannWhitneyUSignificant, + -- * Wilcoxon signed-rank matched-pair test (non-parametric equivalent to the paired t-test) + wilcoxonMatchedPairSignedRank, wilcoxonMatchedPairSignificant, wilcoxonMatchedPairSignificance, wilcoxonMatchedPairCriticalValue, + -- * Wilcoxon rank sum test + wilcoxonRankSums) where + +import Control.Applicative ((<$>)) +import Control.Arrow ((***)) +import Data.Function (on) +import Data.List (findIndex, groupBy, partition, sortBy) +import Data.Ord (comparing) +import qualified Data.Vector.Unboxed as U (length, toList, zipWith) + +import Statistics.Distribution (quantile) +import Statistics.Distribution.Normal (standard) +import Statistics.Math (choose) +import Statistics.Types (Sample) + +-- | The Wilcoxon Rank Sums Test. +-- +-- This test calculates the sum of ranks for the given two samples. The samples +-- are ordered, and assigned ranks (ties are given their average rank), then these +-- ranks are summed for each sample. +-- +-- The return value is (W_1, W_2) where W_1 is the sum of ranks of the first sample +-- and W_2 is the sum of ranks of the second sample. This test is trivially transformed +-- into the Mann-Whitney U test. You will probably want to use 'mannWhitneyU' +-- and the related functions for testing significance, but this function is exposed +-- for completeness. +wilcoxonRankSums :: Sample -> Sample -> (Double, Double) +wilcoxonRankSums xs1 xs2 + = ((sum . map fst) *** (sum . map fst)) . -- sum the ranks per group + partition snd . -- split them back into left and right + concatMap mergeRanks . -- merge the ranks of duplicates + groupBy ((==) `on` (snd . snd)) . -- group duplicate values + zip [1..] . -- give them ranks (duplicates receive different ranks here) + sortBy (comparing snd) $ -- sort by their values + zip (repeat True) (U.toList xs1) ++ zip (repeat False) (U.toList xs2) + -- Tag each sample with an identifier before we merge them + where + mergeRanks :: [(AbsoluteRank, (Bool, Double))] -> [(AbsoluteRank, Bool)] + mergeRanks xs = zip (repeat rank) (map (fst . snd) xs) + where + -- Ranks are merged by assigning them all the average of their ranks: + rank = sum (map fst xs) / fromIntegral (length xs) + +-- | The Mann-Whitney U Test. +-- +-- This is sometimes known as the Mann-Whitney-Wilcoxon U test, and +-- confusingly many sources state that the Mann-Whitney U test is the same as +-- the Wilcoxon's rank sum test (which is provided as 'wilcoxonRankSums'). +-- The Mann-Whitney U is a simple transform of Wilcoxon's rank sum test. +-- +-- Again confusingly, different sources state reversed definitions for U_1 and U_2, +-- so it is worth being explicit about what this function returns. Given two samples, +-- the first, xs_1, of size n_1 and the second, xs_2, of size n_2, this function +-- returns (U_1, U_2) where U_1 = W_1 - (n_1*(n_1+1))\/2 and U_2 = W_2 - (n_2*(n_2+1))\/2, +-- where (W_1, W_2) is the return value of @wilcoxonRankSums xs1 xs2@. +-- +-- Some sources instead state that U_1 and U_2 should be the other way round, often +-- expressing this using U_1' = n_1*n_2 - U_1 (since U_1 + U_2 = n_1*n*2). +-- +-- All of which you probably don't care about if you just feed this into 'mannWhitneyUSignificant'. +mannWhitneyU :: Sample -> Sample -> (Double, Double) +mannWhitneyU xs1 xs2 + = (fst summedRanks - (n1*(n1 + 1))/2 + ,snd summedRanks - (n2*(n2 + 1))/2) + where + n1 = fromIntegral $ U.length xs1 + n2 = fromIntegral $ U.length xs2 + + summedRanks = wilcoxonRankSums xs1 xs2 + +-- | Calculates the critical value of Mann-Whitney U for the given sample +-- sizes and significance level. +-- +-- This function returns the exact calculated value of U for all sample sizes; +-- it does not use the normal approximation at all. Above sample size 20 it is +-- generally recommended to use the normal approximation instead, but this function +-- will calculate the higher critical values if you need them. +-- +-- The algorithm to generate these values is a faster, memoised version of the +-- simple unoptimised generating function given in section 2 of \"The Mann Whitney +-- Wilcoxon Distribution Using Linked Lists\", Cheung and Klotz, Statistica Sinica +-- 7 (1997), . +mannWhitneyUCriticalValue :: (Int, Int) -- ^ The sample size + -> Double -- ^ The p-value (e.g. 0.05) for which you want the critical value. + -> Maybe Int -- ^ The critical value (of U). +mannWhitneyUCriticalValue (m, n) p + | p' <= 1 = Nothing + | m < 1 || n < 1 = Nothing + | otherwise = findIndex (>= p') $ let + firstHalf = map fromIntegral $ take (((m*n)+1)`div`2) $ tail $ alookup !! (m+n-2) !! (min m n - 1) + {- Original: [fromIntegral $ a k (m+n) (min m n) | k <- [1..m*n]] -} + secondHalf + | even (m*n) = reverse firstHalf + | otherwise = tail $ reverse firstHalf + in firstHalf ++ map (mnCn -) secondHalf + where + mnCn = (m+n) `choose` n + p' = mnCn * p + +{- Original function, without memoisation, from Cheung and Klotz: +a :: Int -> Int -> Int -> Int +a u bigN m + | u < 0 = 0 + | u >= (m * smalln) = floor $ fromIntegral bigN `choose` fromIntegral m + | m == 1 || smalln == 1 = u + 1 + | otherwise = a u (bigN - 1) m + + a (u - smalln) (bigN - 1) (m-1) + where smalln = bigN - m +-} + +-- Memoised version of the original a function, above. +-- +-- outer list is indexed by big N - 2 +-- inner list by m (we know m < bigN) +-- innermost list by u +-- +-- So: (alookup ! (bigN - 2) ! m ! u) == a u bigN m +alookup :: [[[Int]]] +alookup = gen 2 [1 : repeat 2] + where + gen bigN predBigNList + = let bigNlist = [ let limit = round $ fromIntegral bigN `choose` fromIntegral m + in [amemoed u m | u <- [0..m*(bigN-m)]] ++ repeat limit + | m <- [1..(bigN-1)]] -- has bigN-1 elements + in bigNlist : gen (bigN+1) bigNlist + where + amemoed :: Int -> Int -> Int + amemoed u m + | m == 1 || smalln == 1 = u + 1 + | otherwise = let (predmList : mList : _) = drop (m-2) predBigNList -- m-2 because starts at 1 + -- We know that predBigNList has bigN - 2 elements + -- (and we know that smalln > 1 therefore bigN > m + 1) + -- So bigN - 2 >= m, i.e. predBigNList must have at least m elements + -- elements, so dropping (m-2) must leave at least 2 + in (mList !! u) + (if u < smalln then 0 else predmList !! (u - smalln)) + where smalln = bigN - m + +-- | Calculates whether the Mann Whitney U test is significant. +-- +-- If both sample sizes are less than or equal to 20, the exact U critical value +-- (as calculated by 'mannWhitneyUCriticalValue') is used. If either sample is +-- larger than 20, the normal approximation is used instead. +-- +-- If you use a one-tailed test, the test indicates whether the first sample is +-- significantly larger than the second. If you want the opposite, simply reverse +-- the order in both the sample size and the (U_1, U_2) pairs. +mannWhitneyUSignificant :: Bool -- ^ Perform one-tailed test (see description above). + -> (Int, Int) -- ^ The sample size from which the (U_1,U_2) values were derived. + -> Double -- ^ The p-value at which to test (e.g. 0.05) + -> (Double, Double) -- ^ The (U_1, U_2) values from 'mannWhitneyU'. + -> Maybe Bool -- ^ Just True if the test is significant, Just + -- False if it is not, and Nothing if the sample + -- was too small to make a decision. +mannWhitneyUSignificant oneTail (in1, in2) p (u1, u2) + | in1 > 20 || in2 > 20 --Use normal approximation +-- = (n1*(n1+1))/2 - u1 - (n1*(n1+n2))/2 +-- = (n1*(n1+1))/2 - (-2*u1 + n1*(n1+n2))/2 +-- = (n1*(n1+1) - 2*u1 + n1*(n1+n2))/2 +-- = (n1*(2*n1 + n2 + 1) - 2*u1)/2 + = let num = (n1*(2*n1 + n2 + 1)) / 2 - u1 + denom = sqrt $ n1*n2*(n1 + n2 + 1) / 12 + z = num / denom + zcrit = quantile standard (1 - if oneTail then p else p/2) + in Just $ (if oneTail then z else abs z) > zcrit + | otherwise = do crit <- fromIntegral <$> mannWhitneyUCriticalValue (in1, in2) p + return $ if oneTail + then u2 <= crit + else min u1 u2 <= crit + where + n1 = fromIntegral in1 + n2 = fromIntegral in2 + +-- | The Wilcoxon matched-pairs signed-rank test. +-- +-- The value returned is the pair (T+, T-). T+ is the sum of positive ranks (the +-- ranks of the differences where the first parameter is higher) whereas T- is +-- the sum of negative ranks (the ranks of the differences where the second parameter is higher). +-- These values mean little by themselves, and should be combined with the 'wilcoxonSignificant' +-- function in this module to get a meaningful result. +-- +-- The samples are zipped together: if one is longer than the other, both are truncated +-- to the the length of the shorter sample. +-- +-- Note that: wilcoxonMatchedPairSignedRank == (\(x, y) -> (y, x)) . flip wilcoxonMatchedPairSignedRank +wilcoxonMatchedPairSignedRank :: Sample -> Sample -> (Double, Double) +wilcoxonMatchedPairSignedRank a b + -- Best to read this function bottom to top: + = (sum *** sum) . -- Sum the positive and negative ranks separately. + partition (> 0) . -- Split the ranks into positive and negative. None of the + -- ranks can be zero. + concatMap mergeRanks . -- Then merge the ranks for any duplicates by taking + -- the average of the ranks, and also make the rank + -- into a signed rank + groupBy ((==) `on` abs . snd) . -- Now group any duplicates together + -- Note: duplicate means same absolute difference + zip [1..] . -- Add a rank (note: at this stage, duplicates will get different ranks) + dropWhile (== 0) . -- Remove any differences that are zero (i.e. ties in the + -- original data). We know they must be at the head of + -- the list because we just sorted it, so dropWhile not filter + sortBy (comparing abs) . -- Sort the differences by absolute difference + U.toList $ -- Convert to a list (could be done later in the pipeline?) + U.zipWith (-) a b -- Work out differences + where + mergeRanks :: [(AbsoluteRank, Double)] -> [SignedRank] + mergeRanks xs = map ((* rank) . signum . snd) xs + -- Note that signum above will always be 1 or -1; any zero differences will + -- have been removed before this function is called. + where + -- Ranks are merged by assigning them all the average of their ranks: + rank = sum (map fst xs) / fromIntegral (length xs) + +type AbsoluteRank = Double +type SignedRank = Double + +-- | The coefficients for x^0, x^1, x^2, etc, in the expression +-- \prod_{r=1}^s (1 + x^r). See the Mitic paper for details. +-- +-- We can define: +-- f(1) = 1 + x +-- f(r) = (1 + x^r)*f(r-1) +-- = f(r-1) + x^r * f(r-1) +-- The effect of multiplying the equation by x^r is to shift +-- all the coefficients by r down the list. +-- +-- This list will be processed lazily from the head. +coefficients :: Int -> [Int] +coefficients 1 = [1, 1] -- 1 + x +coefficients r = let coeffs = coefficients (r-1) + (firstR, rest) = splitAt r coeffs + in firstR ++ add rest coeffs + where + add (x:xs) (y:ys) = x + y : add xs ys + add xs [] = xs + add [] ys = ys + +-- This list will be processed lazily from the head. +summedCoefficients :: Int -> [Double] +summedCoefficients = map fromIntegral . scanl1 (+) . coefficients + +-- | Tests whether a given result from a Wilcoxon signed-rank matched-pairs test +-- is significant at the given level. +-- +-- This function can perform a one-tailed or two-tailed test. If the first +-- parameter to this function is False, the test is performed two-tailed to +-- check if the two samples differ significantly. If the first parameter is +-- True, the check is performed one-tailed to decide whether the first sample +-- (i.e. the first sample you passed to 'wilcoxonMatchedPairSignedRank') is +-- greater than the second sample (i.e. the second sample you passed to +-- 'wilcoxonMatchedPairSignedRank'). If you wish to perform a one-tailed test +-- in the opposite direction, you can either pass the parameters in a different +-- order to 'wilcoxonMatchedPairSignedRank', or simply swap the values in the resulting +-- pair before passing them to this function. +wilcoxonMatchedPairSignificant :: Bool -- ^ Perform one-tailed test (see description above). + -> Int -- ^ The sample size from which the (T+,T-) values were derived. + -> Double -- ^ The p-value at which to test (e.g. 0.05) + -> (Double, Double) -- ^ The (T+, T-) values from 'wilcoxonMatchedPairSignedRank'. + -> Maybe Bool -- ^ Just True if the test is significant, Just + -- False if it is not, and Nothing if the sample + -- was too small to make a decision. +wilcoxonMatchedPairSignificant oneTail sampleSize p (tPlus, tMinus) + -- According to my nearest book (Understanding Research Methods and Statistics + -- by Gary W. Heiman, p590), to check that the first sample is bigger you must + -- use the absolute value of T- for a one-tailed check: + | oneTail = ((abs tMinus <=) . fromIntegral) <$> wilcoxonMatchedPairCriticalValue sampleSize p + -- Otherwise you must use the value of T+ and T- with the smallest absolute value: + | otherwise = ((t <=) . fromIntegral) <$> wilcoxonMatchedPairCriticalValue sampleSize (p/2) + where + t = min (abs tPlus) (abs tMinus) + +-- | Obtains the critical value of T to compare against, given a sample size +-- and a p-value (significance level). Your T value must be less than or +-- equal to the return of this function in order for the test to work out +-- significant. If there is a Nothing return, the sample size is too small to +-- make a decision. +-- +-- 'wilcoxonSignificant' tests the return value of 'wilcoxonMatchedPairSignedRank' +-- for you, so you should use 'wilcoxonSignificant' for determining test results. +-- However, this function is useful, for example, for generating lookup tables +-- for Wilcoxon signed rank critical values. +-- +-- The return values of this function are generated using the method detailed in +-- the paper \"Critical Values for the Wilcoxon Signed Rank Statistic\", Peter +-- Mitic, The Mathematica Journal, volume 6, issue 3, 1996, which can be found +-- here: . +-- According to that paper, the results may differ from other published lookup tables, but +-- (Mitic claims) the values obtained by this function will be the correct ones. +wilcoxonMatchedPairCriticalValue :: Int -- ^ The sample size + -> Double -- ^ The p-value (e.g. 0.05) for which you want the critical value. + -> Maybe Int -- ^ The critical value (of T), or Nothing if + -- the sample is too small to make a decision. +wilcoxonMatchedPairCriticalValue sampleSize p + = case critical of + Just n | n < 0 -> Nothing + | otherwise -> Just n + Nothing -> Just maxBound -- shouldn't happen: beyond end of list + where + m = (2 ** fromIntegral sampleSize) * p + critical = subtract 1 <$> findIndex (> m) (summedCoefficients sampleSize) + +-- | Works out the significance level (p-value) of a T value, given a sample +-- size and a T value from the Wilcoxon signed-rank matched-pairs test. +-- +-- See the notes on 'wilcoxonCriticalValue' for how this is calculated. +wilcoxonMatchedPairSignificance :: Int -- ^ The sample size + -> Double -- ^ The value of T for which you want the significance. + -> Double -- ^^ The significance (p-value). +wilcoxonMatchedPairSignificance sampleSize rank + = (summedCoefficients sampleSize !! floor rank) / 2 ** fromIntegral sampleSize + diff --git a/statistics.cabal b/statistics.cabal --- a/statistics.cabal +++ b/statistics.cabal @@ -1,5 +1,5 @@ name: statistics -version: 0.6.0.2 +version: 0.7.0.0 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful @@ -7,7 +7,7 @@ robustness, and use of good algorithms. Where possible, we provide references to the statistical literature. . - The library's facilities can be divided into three broad categories: + The library's facilities can be divided into four broad categories: . Working with widely used discrete and continuous probability distributions. (There are dozens of exotic distributions in use; we @@ -17,6 +17,8 @@ estimation, bootstrap methods, and autocorrelation analysis. . Random variate generation under several different distributions. + . + Common statistical tests for significant differences between samples. license: BSD3 license-file: LICENSE homepage: http://darcs.serpentine.com/statistics @@ -25,7 +27,7 @@ copyright: 2009, 2010 Bryan O'Sullivan category: Math, Statistics build-type: Simple -cabal-version: >= 1.2 +cabal-version: >= 1.6 extra-source-files: README library @@ -34,6 +36,7 @@ Statistics.Constants Statistics.Distribution Statistics.Distribution.Binomial + Statistics.Distribution.ChiSquared Statistics.Distribution.Gamma Statistics.Distribution.Geometric Statistics.Distribution.Exponential @@ -48,6 +51,7 @@ Statistics.Resampling.Bootstrap Statistics.Sample Statistics.Sample.Powers + Statistics.Test.NonParametric Statistics.Types other-modules: Statistics.Internal @@ -69,3 +73,7 @@ ghc-options: -Wall -funbox-strict-fields if impl(ghc >= 6.8) ghc-options: -fwarn-tabs + +source-repository head + type: mercurial + location: http://bitbucket.org/bos/statistics