statistics 0.6.0.2 → 0.7.0.0
raw patch · 13 files changed
+651/−243 lines, 13 files
Files
- Statistics/Distribution.hs +43/−10
- Statistics/Distribution/Binomial.hs +23/−73
- Statistics/Distribution/ChiSquared.hs +80/−0
- Statistics/Distribution/Exponential.hs +31/−15
- Statistics/Distribution/Gamma.hs +14/−6
- Statistics/Distribution/Geometric.hs +18/−17
- Statistics/Distribution/Hypergeometric.hs +20/−45
- Statistics/Distribution/Normal.hs +24/−24
- Statistics/Distribution/Poisson.hs +24/−33
- Statistics/Math.hs +8/−11
- Statistics/Sample.hs +30/−6
- Statistics/Test/NonParametric.hs +325/−0
- statistics.cabal +11/−3
Statistics/Distribution.hs view
@@ -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 #-}
Statistics/Distribution/Binomial.hs view
@@ -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")
+ Statistics/Distribution/ChiSquared.hs view
@@ -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 #-}
Statistics/Distribution/Exponential.hs view
@@ -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 #-}
Statistics/Distribution/Gamma.hs view
@@ -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
Statistics/Distribution/Geometric.hs view
@@ -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 #-}
Statistics/Distribution/Hypergeometric.hs view
@@ -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 #-}
Statistics/Distribution/Normal.hs view
@@ -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
Statistics/Distribution/Poisson.hs view
@@ -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 #-}
Statistics/Math.hs view
@@ -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))
Statistics/Sample.hs view
@@ -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.
+ Statistics/Test/NonParametric.hs view
@@ -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), <http://www3.stat.sinica.edu.tw/statistica/oldpdf/A7n316.pdf>.+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: <http://www.mathematica-journal.com/issue/v6i3/article/mitic/contents/63mitic.pdf>.+-- 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+
statistics.cabal view
@@ -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