statistics 0.1 → 0.2
raw patch · 19 files changed
+1051/−92 lines, 19 filesdep +mersenne-randomdep ~uvector-algorithms
Dependencies added: mersenne-random
Dependency ranges changed: uvector-algorithms
Files
- Statistics/Autocorrelation.hs +46/−0
- Statistics/Constants.hs +17/−5
- Statistics/Distribution.hs +14/−6
- Statistics/Distribution/Binomial.hs +78/−0
- Statistics/Distribution/Exponential.hs +56/−0
- Statistics/Distribution/Gamma.hs +66/−0
- Statistics/Distribution/Geometric.hs +61/−0
- Statistics/Distribution/Hypergeometric.hs +107/−0
- Statistics/Distribution/Normal.hs +27/−16
- Statistics/Distribution/Poisson.hs +63/−0
- Statistics/Function.hs +6/−6
- Statistics/KernelDensity.hs +4/−0
- Statistics/Math.hs +237/−0
- Statistics/Quantile.hs +42/−34
- Statistics/Resampling.hs +76/−0
- Statistics/Resampling/Bootstrap.hs +95/−0
- Statistics/Sample.hs +19/−19
- Statistics/Types.hs +8/−0
- statistics.cabal +29/−6
+ Statistics/Autocorrelation.hs view
@@ -0,0 +1,46 @@+-- |+-- Module : Statistics.Autocorrelation+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- Functions for computing autocovariance and autocorrelation of a+-- sample.++module Statistics.Autocorrelation+ (+ autocovariance+ , autocorrelation+ ) where++import Data.Array.Vector+import Statistics.Sample (Sample, mean)++-- | Compute the autocovariance of a sample, i.e. the covariance of+-- the sample against a shifted version of itself.+autocovariance :: Sample -> UArr Double+autocovariance a = mapU f . enumFromToU 0 $ l-2+ where+ f k = sumU (zipWithU (*) (takeU (l-k) c) (sliceU c k (l-k)))+ / fromIntegral l+ c = mapU (subtract (mean a)) a+ l = lengthU a++-- | Compute the autocorrelation function of a sample, and the upper+-- and lower bounds of confidence intervals for each element.+--+-- /Note/: The calculation of the 95% confidence interval assumes a+-- stationary Gaussian process.+autocorrelation :: Sample -> (UArr Double, UArr Double, UArr Double)+autocorrelation a = (r, ci (-), ci (+))+ where+ r = mapU (/ headU c) c+ where c = autocovariance a+ dllse = mapU f . scanl1U (+) . mapU square $ r+ where f v = 1.96 * sqrt ((v * 2 + 1) / l)+ l = fromIntegral (lengthU a)+ ci f = consU 1 . tailU . mapU (f (-1/l)) $ dllse+ square x = x * x
Statistics/Constants.hs view
@@ -11,9 +11,11 @@ module Statistics.Constants (- m_huge+ m_epsilon+ , m_huge , m_1_sqrt_2 , m_2_sqrt_pi+ , m_max_exp , m_sqrt_2 , m_sqrt_2_pi ) where@@ -23,22 +25,32 @@ m_huge = 1.797693e308 {-# INLINE m_huge #-} +-- | The largest 'Int' /x/ such that 2**(/x/-1) is approximately+-- representable as a 'Double'.+m_max_exp :: Int+m_max_exp = 1024+ -- | @sqrt 2@ m_sqrt_2 :: Double-m_sqrt_2 = 1.414213562373095145474621858739+m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766 {-# INLINE m_sqrt_2 #-} -- | @sqrt (2 * pi)@ m_sqrt_2_pi :: Double-m_sqrt_2_pi = 2.506628274631000241612355239340+m_sqrt_2_pi = 2.5066282746310005024157652848110452530069867406099383166299 {-# INLINE m_sqrt_2_pi #-} -- | @2 / sqrt pi@ m_2_sqrt_pi :: Double-m_2_sqrt_pi = 1.128379167095512558560699289956+m_2_sqrt_pi = 1.1283791670955125738961589031215451716881012586579977136881 {-# INLINE m_2_sqrt_pi #-} -- | @1 / sqrt 2@ m_1_sqrt_2 :: Double-m_1_sqrt_2 = 0.707106781186547461715008466854+m_1_sqrt_2 = 0.7071067811865475244008443621048490392848359376884740365883 {-# INLINE m_1_sqrt_2 #-}++-- | The smallest 'Double' larger than 1.+m_epsilon :: Double+m_epsilon = encodeFloat (signif+1) expo - 1.0+ where (signif,expo) = decodeFloat (1.0::Double)
Statistics/Distribution.hs view
@@ -13,31 +13,39 @@ module Statistics.Distribution ( Distribution(..)+ , Mean(..)+ , Variance(..) , findRoot ) where -- | The interface shared by all probability distributions. class Distribution d where -- | Probability density function. The probability that a- -- stochastic variable @x@ has the value @X@, i.e. @P(x=X)@.+ -- stochastic variable /x/ has the value /X/, i.e. P(/x/=/X/). probability :: d -> Double -> Double -- | Cumulative distribution function. The probability that a- -- stochastic variable @x@ is less than @X@, i.e. @P(x<X)@.+ -- stochastic variable /x/ is less 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)@.+ -- /X/ for which P(/x/</X/). inverse :: d -> Double -> Double --- | Approximate the value of @X@ for which @P(x>X) == p@.+class Distribution d => Mean d where+ mean :: d -> Double++class Mean d => Variance d where+ variance :: d -> Double++-- | Approximate the value of /X/ for which P(/x/>/X/)=/p/. -- -- This method uses a combination of Newton-Raphson iteration and -- 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@.+-- distribution reaches the value /p/. findRoot :: Distribution d => d- -> Double -- ^ Probability @p@+ -> Double -- ^ Probability /p/ -> Double -- ^ Initial guess -> Double -- ^ Lower bound on interval -> Double -- ^ Upper bound on interval
+ Statistics/Distribution/Binomial.hs view
@@ -0,0 +1,78 @@+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Statistics.Distribution.Binomial+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The binomial distribution. This is the discrete probability+-- distribution of the number of successes in a sequence of /n/+-- independent yes\/no experiments, each of which yields success with+-- probability /p/.++module Statistics.Distribution.Binomial+ (+ BinomialDistribution+ -- * Constructors+ , binomial+ -- * Accessors+ , bdTrials+ , bdProbability+ ) where++import Control.Exception (assert)+import Data.Array.Vector+import Data.Typeable (Typeable)+import qualified Statistics.Distribution as D+import Statistics.Math (choose)++-- | The binomial distribution.+data BinomialDistribution = BD {+ bdTrials :: {-# UNPACK #-} !Int+ -- ^ Number of trials.+ , bdProbability :: {-# UNPACK #-} !Double+ -- ^ Probability.+ } deriving (Eq, Read, Show, Typeable)++instance D.Distribution BinomialDistribution where+ probability = probability+ cumulative = cumulative+ inverse = inverse++instance D.Variance BinomialDistribution where+ variance = variance++instance D.Mean BinomialDistribution where+ mean = mean++probability :: BinomialDistribution -> Double -> Double+probability (BD n p) x =+ fromIntegral (n `choose` floor x) * p ** x * (1-p) ** (fromIntegral n-x)++cumulative :: BinomialDistribution -> Double -> Double+cumulative d =+ sumU . mapU (probability d . fromIntegral) . enumFromToU (0::Int) . floor++inverse :: BinomialDistribution -> Double -> Double+inverse d@(BD n _p) p = D.findRoot d p (n'/2) 0 n'+ where n' = fromIntegral n++mean :: BinomialDistribution -> Double+mean (BD n p) = fromIntegral n * p+{-# INLINE mean #-}++variance :: BinomialDistribution -> Double+variance (BD n p) = fromIntegral n * p * (1 - p)+{-# INLINE variance #-}++binomial :: Int -- ^ Number of trials.+ -> Double -- ^ Probability.+ -> BinomialDistribution+binomial n p =+ assert (n > 0) .+ assert (p > 0 && p < 1) $+ BD n p+{-# INLINE binomial #-}
+ Statistics/Distribution/Exponential.hs view
@@ -0,0 +1,56 @@+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Statistics.Distribution.Exponential+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The exponential distribution. This is the discrete probability+-- distribution of the number of successes in a sequence of /n/+-- independent yes\/no experiments, each of which yields success with+-- probability /p/.++module Statistics.Distribution.Exponential+ (+ ExponentialDistribution+ -- * Constructors+ , fromLambda+ , fromSample+ ) where++import Data.Typeable (Typeable)+import qualified Statistics.Distribution as D+import qualified Statistics.Sample as S+import Statistics.Types (Sample)++newtype ExponentialDistribution = ED {+ edLambda :: Double+ } deriving (Eq, Read, Show, Typeable)++instance D.Distribution ExponentialDistribution where+ probability (ED l) x = l * exp (-l * x)+ {-# INLINE probability #-}+ cumulative (ED l) x = 1 - exp (-l * x)+ {-# INLINE cumulative #-}+ inverse (ED l) p = -log (1 - p) / l+ {-# INLINE inverse #-}++instance D.Variance ExponentialDistribution where+ variance (ED l) = l * l+ {-# INLINE variance #-}++instance D.Mean ExponentialDistribution where+ mean = edLambda+ {-# INLINE mean #-}++fromLambda :: Double -- ^ λ (scale) parameter.+ -> ExponentialDistribution+fromLambda = ED+{-# INLINE fromLambda #-}++fromSample :: Sample -> ExponentialDistribution+fromSample = ED . S.mean+{-# INLINE fromSample #-}
+ Statistics/Distribution/Gamma.hs view
@@ -0,0 +1,66 @@+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Statistics.Distribution.Gamma+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The gamma distribution. This is a continuous probability+-- distribution with two parameters, /k/ and ϑ. If /k/ is+-- integral, the distribution represents the sum of /k/ independent+-- exponentially distributed random variables, each of which has a+-- mean of ϑ.++module Statistics.Distribution.Gamma+ (+ GammaDistribution+ -- * Constructors+ --, fromParams+ --, fromSample+ --, standard+ -- * Accessors+ , gdShape+ , gdScale+ ) where++import Data.Typeable (Typeable)+import Statistics.Constants (m_huge)+import Statistics.Math (incompleteGamma, logGamma)+import qualified Statistics.Distribution as D++-- | The gamma distribution.+data GammaDistribution = GD {+ gdShape :: {-# UNPACK #-} !Double -- ^ Shape parameter, /k/.+ , gdScale :: {-# UNPACK #-} !Double -- ^ Scale parameter, ϑ.+ } deriving (Eq, Read, Show, Typeable)++instance D.Distribution GammaDistribution where+ probability = probability+ cumulative = cumulative+ inverse = inverse++instance D.Variance GammaDistribution where+ variance (GD a l) = a / (l * l)+ {-# INLINE variance #-}++instance D.Mean GammaDistribution where+ mean (GD a l) = a / l+ {-# INLINE mean #-}++probability :: GammaDistribution -> Double -> Double+probability (GD a l) x = x ** (a-1) * exp (-x/l) / (exp (logGamma a) * l ** a)+{-# INLINE probability #-}++cumulative :: GammaDistribution -> Double -> Double+cumulative (GD a l) x = incompleteGamma a (x/l) / exp (logGamma a)+{-# INLINE cumulative #-}++inverse :: GammaDistribution -> Double -> Double+inverse d p+ | p == 0 = -1/0+ | p == 1 = 1/0+ | otherwise = D.findRoot d p (gdShape d) 0 m_huge+{-# INLINE inverse #-}
+ Statistics/Distribution/Geometric.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Statistics.Distribution.Geometric+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The Geometric distribution. This is the discrete probability+-- distribution of a number of events occurring in a fixed interval if+-- these events occur with a known average rate, and occur+-- independently from each other within that interval.++module Statistics.Distribution.Geometric+ (+ GeometricDistribution+ -- * Constructors+ , fromSuccess+ -- ** Accessors+ , pdSuccess+ ) where++import Control.Exception (assert)+import Data.Typeable (Typeable)+import qualified Statistics.Distribution as D++newtype GeometricDistribution = GD {+ pdSuccess :: Double+ } deriving (Eq, Read, Show, Typeable)++instance D.Distribution GeometricDistribution where+ probability = probability+ cumulative = cumulative+ inverse = inverse++instance D.Variance GeometricDistribution where+ variance (GD s) = (1 - s) / (s * s)+ {-# INLINE variance #-}++instance D.Mean GeometricDistribution where+ mean (GD s) = 1 / s+ {-# INLINE mean #-}++fromSuccess :: Double -> GeometricDistribution+fromSuccess x = assert (x >= 0 && x <= 1)+ GD x+{-# INLINE fromSuccess #-}++probability :: GeometricDistribution -> Double -> Double+probability (GD s) x = s * (1-s) ** (x-1)+{-# INLINE probability #-}++cumulative :: GeometricDistribution -> Double -> Double+cumulative (GD s) x = 1 - (1-s) ** x+{-# INLINE cumulative #-}++inverse :: GeometricDistribution -> Double -> Double+inverse (GD s) p = log (1 - p) / log (1 - s)+{-# INLINE inverse #-}
+ Statistics/Distribution/Hypergeometric.hs view
@@ -0,0 +1,107 @@+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Statistics.Distribution.Hypergeometric+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The Hypergeometric distribution. This is the discrete probability+-- distribution that measures the probability of /k/ successes in /l/+-- trials, without replacement, from a finite population.+--+-- The parameters of the distribution describe /k/ elements chosen+-- from a population of /l/, with /m/ elements of one type, and+-- /l/-/m/ of the other (all are positive integers).++module Statistics.Distribution.Hypergeometric+ (+ HypergeometricDistribution+ -- * Constructors+ , fromParams+ -- ** Accessors+ , hdM+ , hdL+ , hdK+ ) where++import Control.Exception (assert)+import Data.Array.Vector+import Data.Typeable (Typeable)+import Statistics.Math (choose, logFactorial)+import Statistics.Constants (m_max_exp)+import qualified Statistics.Distribution as D++data HypergeometricDistribution = HD {+ hdM :: {-# UNPACK #-} !Int+ , hdL :: {-# UNPACK #-} !Int+ , hdK :: {-# UNPACK #-} !Int+ } deriving (Eq, Read, Show, Typeable)++instance D.Distribution HypergeometricDistribution where+ probability = probability+ cumulative = cumulative+ inverse = inverse++instance D.Variance HypergeometricDistribution where+ variance = variance++instance D.Mean HypergeometricDistribution where+ mean = mean++variance :: HypergeometricDistribution -> Double+variance (HD m l k) = (k' * ml) * (1 - ml) * (l' - k') / (l' - 1)+ where m' = fromIntegral m+ l' = fromIntegral l+ k' = fromIntegral k+ ml = m' / l'+{-# INLINE variance #-}++mean :: HypergeometricDistribution -> Double+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) .+ assert (l > 0) .+ assert (k > 0 && k <= l) $+ HD m l k+{-# INLINE fromParams #-}++probability :: HypergeometricDistribution -> Double -> Double+probability (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 = fromIntegral (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 probability #-}++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 = sumU . mapU (probability d . fromIntegral) . enumFromToU imin . floor $ x+{-# INLINE cumulative #-}++inverse :: HypergeometricDistribution -> Double -> Double+inverse = error "Statistics.Distribution.Hypergeometric.inverse: not yet implemented"+{-# INLINE inverse #-}
Statistics/Distribution/Normal.hs view
@@ -1,5 +1,6 @@+{-# LANGUAGE DeriveDataTypeable #-} -- |--- Module : Statistics.Normal+-- Module : Statistics.Distribution.Normal -- Copyright : (c) 2009 Bryan O'Sullivan -- License : BSD3 --@@ -7,11 +8,13 @@ -- Stability : experimental -- Portability : portable ----- The normal distribution.+-- The normal distribution. This is a continuous probability+-- distribution that describes data that cluster around a mean. module Statistics.Distribution.Normal ( NormalDistribution+ -- * Constructors , fromParams , fromSample , standard@@ -19,50 +22,58 @@ import Control.Exception (assert) import Data.Number.Erf (erfc)+import Data.Typeable (Typeable) import Statistics.Constants (m_huge, m_sqrt_2, m_sqrt_2_pi) import Statistics.Types (Sample) import qualified Statistics.Distribution as D import qualified Statistics.Sample as S -data NormalDistribution = NormalDistribution {+-- | The normal distribution.+data NormalDistribution = ND { mean :: {-# UNPACK #-} !Double , variance :: {-# UNPACK #-} !Double- , pdfDenom :: {-# UNPACK #-} !Double- , cdfDenom :: {-# UNPACK #-} !Double- } deriving (Eq, Ord, Read, Show)+ , ndPdfDenom :: {-# UNPACK #-} !Double+ , ndCdfDenom :: {-# UNPACK #-} !Double+ } deriving (Eq, Read, Show, Typeable) instance D.Distribution NormalDistribution where probability = probability cumulative = cumulative inverse = inverse +instance D.Variance NormalDistribution where+ variance = variance++instance D.Mean NormalDistribution where+ mean = mean+ standard :: NormalDistribution-standard = NormalDistribution {+standard = ND { mean = 0.0 , variance = 1.0- , cdfDenom = m_sqrt_2- , pdfDenom = m_sqrt_2_pi+ , ndPdfDenom = m_sqrt_2_pi+ , ndCdfDenom = m_sqrt_2 } fromParams :: Double -> Double -> NormalDistribution-fromParams m v = assert (v > 0) $- NormalDistribution {+fromParams m v = assert (v > 0)+ ND { mean = m , variance = v- , cdfDenom = m_sqrt_2 * sv- , pdfDenom = m_sqrt_2_pi * sv+ , ndPdfDenom = m_sqrt_2_pi * sv+ , ndCdfDenom = m_sqrt_2 * sv } where sv = sqrt v- + fromSample :: Sample -> NormalDistribution fromSample a = fromParams (S.mean a) (S.variance a) probability :: NormalDistribution -> Double -> Double-probability d x = exp (-xm * xm / (2 * variance d)) / pdfDenom d+probability 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) / cdfDenom d) / 2+cumulative d x = erfc (-(x-mean d) / ndCdfDenom d) / 2 inverse :: NormalDistribution -> Double -> Double inverse d p
+ Statistics/Distribution/Poisson.hs view
@@ -0,0 +1,63 @@+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Statistics.Distribution.Poisson+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The Poisson distribution. This is the discrete probability+-- distribution of a number of events occurring in a fixed interval if+-- these events occur with a known average rate, and occur+-- independently from each other within that interval.++module Statistics.Distribution.Poisson+ (+ PoissonDistribution+ -- * Constructors+ , fromLambda+ -- , fromSample+ ) where++import Data.Array.Vector+import Data.Typeable (Typeable)+import qualified Statistics.Distribution as D+import Statistics.Constants (m_huge)+import Statistics.Math (logGamma)++newtype PoissonDistribution = PD {+ pdLambda :: Double+ } deriving (Eq, Read, Show, Typeable)++instance D.Distribution PoissonDistribution where+ probability = probability+ cumulative = cumulative+ inverse = inverse++instance D.Variance PoissonDistribution where+ variance = pdLambda+ {-# INLINE variance #-}++instance D.Mean PoissonDistribution where+ mean = pdLambda+ {-# INLINE mean #-}++fromLambda :: Double -> PoissonDistribution+fromLambda = PD+{-# INLINE fromLambda #-}++probability :: PoissonDistribution -> Double -> Double+probability (PD l) x = exp (x * log l - l - logGamma x)+{-# INLINE probability #-}++cumulative :: PoissonDistribution -> Double -> Double+cumulative d = sumU . mapU (probability d . fromIntegral) .+ enumFromToU (0::Int) . floor+{-# INLINE cumulative #-}++inverse :: PoissonDistribution -> Double -> Double+inverse d p = fromIntegral . r $ D.findRoot d p (pdLambda d) 0 m_huge+ where r = round :: Double -> Int+{-# INLINE inverse #-}
Statistics/Function.hs view
@@ -1,6 +1,6 @@ {-# LANGUAGE TypeOperators #-} -- |--- Module : Statistics.Quantile+-- Module : Statistics.Function -- Copyright : (c) 2009 Bryan O'Sullivan -- License : BSD3 --@@ -8,7 +8,7 @@ -- Stability : experimental -- Portability : portable ----- Functions for computing quantiles.+-- Useful functions. module Statistics.Function (@@ -17,19 +17,19 @@ , partialSort ) where -import Data.Array.Vector.Algorithms.Immutable (apply)+import Data.Array.Vector.Algorithms.Combinators (apply) import Data.Array.Vector ((:*:)(..), UA, UArr, foldlU) import qualified Data.Array.Vector.Algorithms.Intro as I --- | Sort.+-- | Sort an array. sort :: (UA e, Ord e) => UArr e -> UArr e sort = apply I.sort {-# INLINE sort #-} --- | Partially sort, such that the least @k@ elements will be+-- | Partially sort an array, such that the least /k/ elements will be -- at the front. partialSort :: (UA e, Ord e) =>- Int -- ^ The number @k@ of least elements+ Int -- ^ The number /k/ of least elements. -> UArr e -> UArr e partialSort k = apply (\a -> I.partialSort a k)
Statistics/KernelDensity.hs view
@@ -87,9 +87,13 @@ n' = n - 1 -- | The convolution kernel. Its parameters are as follows:+-- -- * Scaling factor, 1\//nh/+-- -- * Bandwidth, /h/+-- -- * A point at which to sample the input, /p/+-- -- * One sample value, /v/ type Kernel = Double -> Double
+ Statistics/Math.hs view
@@ -0,0 +1,237 @@+{-# LANGUAGE BangPatterns #-}+-- |+-- Module : Statistics.Math+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- Mathematical functions for statistics.++module Statistics.Math+ (+ -- * Functions+ chebyshev+ , choose+ -- ** Factorial functions+ , factorial+ , logFactorial+ -- ** Gamma functions+ , incompleteGamma+ , logGamma+ , logGammaL+ -- * References+ -- $references+ ) where++import Data.Array.Vector+import Data.Word (Word64)+import Statistics.Constants (m_sqrt_2_pi)+import Statistics.Distribution (cumulative)+import Statistics.Distribution.Normal (standard)++data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double++-- | Evaluate a series of Chebyshev polynomials. Uses Clenshaw's+-- algorithm.+chebyshev :: Double -- ^ Parameter of each function.+ -> UArr Double -- ^ Coefficients of each polynomial+ -- term, in increasing order.+ -> Double+chebyshev x a = fini . foldlU step (C 0 0 0) .+ enumFromThenToU (lengthU a - 1) (-1) $ 0+ where step (C u v w) k = C (x2 * v - w + indexU a k) u v+ fini (C u _ w) = (u - w) / 2+ x2 = x * 2++-- | The binomial coefficient.+--+-- > 7 `choose` 3 == 35+choose :: Int -> Int -> Int+n `choose` k+ | k > n = 0+ | otherwise = ceiling . foldlU go 1 . enumFromToU 1 $ k'+ where go a i = a * (nk + j) / j+ where j = fromIntegral i :: Double+ k' | k > n `div` 2 = n - k+ | otherwise = k+ nk = fromIntegral (n - k')+{-# INLINE choose #-}++data F = F {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64++-- | Compute the factorial function /n/!. Returns ∞ if the+-- input is above 170 (above which the result cannot be represented by+-- a 64-bit 'Double').+factorial :: Int -> Double+factorial n+ | n < 0 = error "Statistics.Math.factorial: negative input"+ | n <= 1 = 0+ | n <= 14 = fini . foldlU goLong (F 1 1) $ ns+ | otherwise = foldlU goDouble 1 $ ns+ where goDouble t k = t * fromIntegral k+ goLong (F z x) _ = F (z * x') x'+ where x' = x + 1+ fini (F z _) = fromIntegral z+ ns = enumFromToU 2 n+{-# INLINE factorial #-}++-- | Compute the natural logarithm of the factorial function. Gives+-- 16 decimal digits of precision.+logFactorial :: Int -> Double+logFactorial n+ | n <= 14 = log (factorial n)+ | otherwise = (x - 0.5) * log x - x + 9.1893853320467e-1 + z / x+ where x = fromIntegral (n + 1)+ y = 1 / (x * x)+ z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y -+ 2.7777777777778e-3) * y + 8.3333333333333e-2+{-# INLINE logFactorial #-}++-- | Compute the incomplete gamma integral function γ(/s/,/x/).+-- Uses Algorithm AS 239 by Shea.+incompleteGamma :: Double -- ^ /s/+ -> Double -- ^ /x/+ -> Double+incompleteGamma x p+ | x < 0 || p <= 0 = 1/0+ | x == 0 = 0+ | p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))+ | x >= 1e8 = 0+ | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1)+ g = a + log (pearson p 1 1)+ in if g > limit then exp g else 0+ | otherwise = let g = p * log x - x - logGamma p + log cf+ in if g > limit then 1 - exp g else 1+ where+ norm = cumulative standard+ pearson !a !c !g+ | c' <= tolerance = g'+ | otherwise = pearson a' c' g'+ where a' = a + 1+ c' = c * x / a'+ g' = g + c'+ cf = let a = 1 - p+ b = a + x + 1+ p3 = x + 1+ p4 = x * b+ in contFrac a b 0 1 x p3 p4 (p3/p4)+ contFrac !a !b !c !p1 !p2 !p3 !p4 !g+ | abs (g - rn) <= min tolerance (tolerance * rn) = g+ | otherwise = contFrac a' b' c' (f p3) (f p4) (f p5) (f p6) rn+ where a' = a + 1+ b' = b + 2+ c' = c + 1+ an = a' * c'+ p5 = b' * p3 - an * p1+ p6 = b' * p4 - an * p2+ rn = p5 / p6+ f n | abs p5 > overflow = n / overflow+ | otherwise = n+ limit = -88+ tolerance = 1e-14+ overflow = 1e37++-- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html++-- | Compute the logarithm of the gamma function Γ(/x/). Uses+-- Algorithm AS 245 by Macleod.+--+-- Gives an accuracy of 10–12 significant decimal digits, except+-- for small regions around /x/ = 1 and /x/ = 2, where the function+-- goes to zero. For greater accuracy, use 'logGammaL'.+--+-- Returns ∞ if the input is outside of the range (0 < /x/+-- ≤ 1e305).+logGamma :: Double -> Double+logGamma x+ | x <= 0 = 1/0+ | x < 1.5 = a + c *+ ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) /+ ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)+ | x < 4 = (x - 2) *+ ((((r2_4 * x + r2_3) * x + r2_2) * x + r2_1) * x + r2_0) /+ ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5)+ | x < 12 = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) /+ ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5)+ | x > 5.1e5 = k+ | otherwise = k + x1 *+ ((r4_2 * x2 + r4_1) * x2 + r4_0) /+ ((x2 + r4_4) * x2 + r4_3)+ where+ a :*: b :*: c+ | x < 0.5 = -y :*: x + 1 :*: x+ | otherwise = 0 :*: x :*: x - 1++ y = log x+ k = x * (y-1) - 0.5 * y + alr2pi+ alr2pi = 0.918938533204673++ x1 = 1 / x+ x2 = x1 * x1++ r1_0 = -2.66685511495; r1_1 = -24.4387534237; r1_2 = -21.9698958928+ r1_3 = 11.1667541262; r1_4 = 3.13060547623; r1_5 = 0.607771387771+ r1_6 = 11.9400905721; r1_7 = 31.4690115749; r1_8 = 15.2346874070++ r2_0 = -78.3359299449; r2_1 = -142.046296688; r2_2 = 137.519416416+ r2_3 = 78.6994924154; r2_4 = 4.16438922228; r2_5 = 47.0668766060+ r2_6 = 313.399215894; r2_7 = 263.505074721; r2_8 = 43.3400022514++ r3_0 = -2.12159572323; r3_1 = 2.30661510616; r3_2 = 2.74647644705+ r3_3 = -4.02621119975; r3_4 = -2.29660729780; r3_5 = -1.16328495004+ r3_6 = -1.46025937511; r3_7 = -2.42357409629; r3_8 = -5.70691009324++ r4_0 = 0.279195317918525; r4_1 = 0.4917317610505968;+ r4_2 = 0.0692910599291889; r4_3 = 3.350343815022304+ r4_4 = 6.012459259764103++data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double++-- | Compute the logarithm of the gamma function, Γ(/x/). Uses a+-- Lanczos approximation.+--+-- This function is slower than 'logGamma', but gives 14 or more+-- significant decimal digits of accuracy, except around /x/ = 1 and+-- /x/ = 2, where the function goes to zero.+--+-- Returns ∞ if the input is outside of the range (0 < /x/+-- ≤ 1e305).+logGammaL :: Double -> Double+logGammaL x+ | x <= 0 = 1/0+ | otherwise = fini . foldlU go (L 0 (x+7)) $ a+ where fini (L l _) = log (l+a0) + log m_sqrt_2_pi - x65 + (x-0.5) * log x65+ go (L l t) k = L (l + k / t) (t-1)+ x65 = x + 6.5+ a0 = 0.9999999999995183+ a = toU [ 0.1659470187408462e-06+ , 0.9934937113930748e-05+ , -0.1385710331296526+ , 12.50734324009056+ , -176.6150291498386+ , 771.3234287757674+ , -1259.139216722289+ , 676.5203681218835+ ]++-- $references+--+-- * Clenshaw, C.W. (1962) Chebyshev series for mathematical+-- functions. /National Physical Laboratory Mathematical Tables 5/,+-- Her Majesty's Stationery Office, London.+--+-- * Lanczos, C. (1964) A precision approximation of the gamma+-- function. /SIAM Journal on Numerical Analysis B/+-- 1:86–96. <http://www.jstor.org/stable/2949767>+--+-- * Macleod, A.J. (1989) Algorithm AS 245: A robust and reliable+-- algorithm for the logarithm of the gamma function.+-- /Journal of the Royal Statistical Society, Series C (Applied Statistics)/+-- 38(2):397–402. <http://www.jstor.org/stable/2348078>+--+-- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete+-- gamma integral. /Applied Statistics/+-- 37(3):466–473. <http://www.jstor.org/stable/2347328>
Statistics/Quantile.hs view
@@ -8,15 +8,20 @@ -- Stability : experimental -- Portability : portable ----- Functions for approximating quantiles.+-- Functions for approximating quantiles, i.e. points taken at regular+-- intervals from the cumulative distribution function of a random+-- variable.+--+-- The number of quantiles is described below by the variable /q/, so+-- with /q/=4, a 4-quantile (also known as a /quartile/) has 4+-- intervals, and contains 5 points. The parameter /k/ describes the+-- desired point, where 0 ≤ /k/ ≤ /q/. module Statistics.Quantile (- -- * Types- ContParam(..)- -- * Quantile estimation functions- , weightedAvg+ weightedAvg+ , ContParam(..) , continuousBy -- * Parameters for the continuous sample method@@ -33,14 +38,15 @@ import Control.Exception (assert) import Data.Array.Vector (allU, indexU, lengthU)+import Statistics.Constants (m_epsilon) import Statistics.Function (partialSort) import Statistics.Types (Sample) --- | Use the weighted average method to estimate the @k@th--- @q@-quantile of a sample.-weightedAvg :: Int -- ^ @k@, the desired quantile- -> Int -- ^ @q@, the number of quantiles- -> Sample -- ^ @x@, the sample data+-- | Estimate the /k/th /q/-quantile of a sample, using the weighted+-- average method.+weightedAvg :: Int -- ^ /k/, the desired quantile.+ -> Int -- ^ /q/, the number of quantiles.+ -> Sample -- ^ /x/, the sample data. -> Double weightedAvg k q x = assert (q >= 2) .@@ -57,15 +63,17 @@ sx = partialSort (j+2) x {-# INLINE weightedAvg #-} --- | Parameters @a@ and @b@ to the 'quantileBy' function.+-- | Parameters /a/ and /b/ to the 'continuousBy' function. data ContParam = ContParam {-# UNPACK #-} !Double {-# UNPACK #-} !Double --- | Using the continuous sample method with the given parameters,--- estimate the @k@th @q@-quantile of a sample @x@.-continuousBy :: ContParam -- ^ Parameters @a@ and @b@- -> Int -- ^ @k@, the desired quantile- -> Int -- ^ @q@, the number of quantiles- -> Sample -- ^ @x@, the sample data+-- | Estimate the /k/th /q/-quantile of a sample /x/, using the+-- continuous sample method with the given parameters. This is the+-- method used by most statistical software, such as R, Mathematica,+-- SPSS, and S.+continuousBy :: ContParam -- ^ Parameters /a/ and /b/.+ -> Int -- ^ /k/, the desired quantile.+ -> Int -- ^ /q/, the number of quantiles.+ -> Sample -- ^ /x/, the sample data. -> Double continuousBy (ContParam a b) k q x = assert (q >= 2) .@@ -80,50 +88,51 @@ h | abs r < eps = 0 | otherwise = r where r = t - fromIntegral j- eps = 8.881784e-16+ eps = m_epsilon * 4 n = lengthU x- item m = indexU sx $ bracket m+ item = indexU sx . bracket sx = partialSort (bracket j + 1) x bracket m = min (max m 0) (n - 1) {-# INLINE continuousBy #-} --- | California Department of Public Works definition, @a=0,b=1@.--- Gives a linear interpolation of the empirical CDF.--- This corresponds to method 4 in R and Mathematica.+-- | California Department of Public Works definition, /a/=0, /b/=1.+-- Gives a linear interpolation of the empirical CDF. This+-- corresponds to method 4 in R and Mathematica. cadpw :: ContParam cadpw = ContParam 0 1 {-# INLINE cadpw #-} --- | Hazen's definition, @a=0.5,b=0.5@. This is claimed to be popular--- among hydrologists. This corresponds to method 5 in R and+-- | Hazen's definition, /a/=0.5, /b/=0.5. This is claimed to be+-- popular among hydrologists. This corresponds to method 5 in R and -- Mathematica. hazen :: ContParam hazen = ContParam 0.5 0.5 {-# INLINE hazen #-} --- | SPSS definition, @a=0,b=0@, also known as Weibull's definition.--- This corresponds to method 6 in R and Mathematica.+-- | Definition used by the SPSS statistics application, with /a/=0,+-- /b/=0 (also known as Weibull's definition). This corresponds to+-- method 6 in R and Mathematica. spss :: ContParam spss = ContParam 0 0 {-# INLINE spss #-} --- | S definition, @a=1,b=1@. The interpolation points divide the--- sample range into @n-1@ intervals. This corresponds to method 7 in--- R and Mathematica.+-- | Definition used by the S statistics application, with /a/=1,+-- /b/=1. The interpolation points divide the sample range into @n-1@+-- intervals. This corresponds to method 7 in R and Mathematica. s :: ContParam s = ContParam 1 1 {-# INLINE s #-} --- | Median unbiased definition, @a=1/3,b=1/3@. The resulting quantile--- estimates are approximately median unbiased regardless of the--- distribution of @x@. This corresponds to method 8 in R and+-- | Median unbiased definition, /a/=1\/3, /b/=1\/3. The resulting+-- quantile estimates are approximately median unbiased regardless of+-- the distribution of /x/. This corresponds to method 8 in R and -- Mathematica. medianUnbiased :: ContParam medianUnbiased = ContParam third third where third = 1/3 {-# INLINE medianUnbiased #-} --- | Normal unbiased definition, @a=3/8,b=3/8@. An approximately+-- | Normal unbiased definition, /a/=3\/8, /b/=3\/8. An approximately -- unbiased estimate if the empirical distribution approximates the -- normal distribution. This corresponds to method 9 in R and -- Mathematica.@@ -140,4 +149,3 @@ -- * Hyndman, R.J.; Fan, Y. (1996) Sample quantiles in statistical -- packages. /American Statistician/ -- 50(4):361–365. <http://www.jstor.org/stable/2684934>-
+ Statistics/Resampling.hs view
@@ -0,0 +1,76 @@+-- |+-- Module : Statistics.Resampling+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- Resampling statistics.++module Statistics.Resampling+ (+ Resample(..)+ , jackknife+ , resample+ ) where++import Control.Exception (assert)+import Control.Monad (forM_)+import Control.Monad.ST (unsafeSTToIO)+import Data.Array.Vector+import Data.Array.Vector.Algorithms.Intro (sort)+import Statistics.Types (Estimator, Sample)+import System.Random.Mersenne (MTGen, random)++-- | A resample drawn randomly, with replacement, from a set of data+-- points. Distinct from a normal array to make it harder for your+-- humble author's brain to go wrong.+newtype Resample = Resample {+ fromResample :: UArr Double+ } deriving (Eq, Show)++-- | Resample a data set repeatedly, with replacement, computing each+-- estimate over the resampled data.+resample :: MTGen -> [Estimator] -> Int -> Sample -> IO [Resample]+resample gen ests numResamples samples = do+ results <- unsafeSTToIO . mapM (const (newMU numResamples)) $ ests+ loop 0 (zip ests results)+ unsafeSTToIO $ do+ mapM_ sort results+ mapM (fmap Resample . unsafeFreezeAllMU) results+ where+ loop k ers | k >= numResamples = return ()+ | otherwise = do+ re <- createU n $ \_ -> do+ r <- random gen+ return (indexU samples (abs r `mod` n))+ unsafeSTToIO . forM_ ers $ \(est,arr) ->+ writeMU arr k . est $ re+ loop (k+1) ers+ n = lengthU samples++-- | Create an array, using the given action to populate each element.+createU :: (UA e) => Int -> (Int -> IO e) -> IO (UArr e)+createU size itemAt = assert (size >= 0) $+ unsafeSTToIO (newMU size) >>= loop 0+ where+ loop k arr | k >= size = unsafeSTToIO (unsafeFreezeAllMU arr)+ | otherwise = do+ r <- itemAt k+ unsafeSTToIO (writeMU arr k r)+ loop (k+1) arr++-- | Compute a statistical estimate repeatedly over a sample, each+-- time omitting a successive element.+jackknife :: Estimator -> Sample -> UArr Double+jackknife est sample = mapU f . enumFromToU 0 . subtract 1 . lengthU $ sample+ where f i = est (dropAt i sample)+{-# INLINE jackknife #-}++-- | Drop the /k/th element of a vector.+dropAt :: UA e => Int -> UArr e -> UArr e+dropAt n = mapU sndT . filterU notN . indexedU+ where notN (i :*: _) = i /= n+ sndT (_ :*: k) = k
+ Statistics/Resampling/Bootstrap.hs view
@@ -0,0 +1,95 @@+-- |+-- Module : Statistics.Resampling.Bootstrap+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- The bootstrap method for statistical inference.++module Statistics.Resampling.Bootstrap+ (+ Estimate(..)+ , bootstrapBCA+ -- * References+ -- $references+ ) where++import Control.Exception (assert)+import Data.Array.Vector (foldlU, filterU, indexU, lengthU)+import Statistics.Distribution.Normal+import Statistics.Distribution (cumulative, inverse)+import Statistics.Resampling (Resample(..), jackknife)+import Statistics.Sample (mean)+import Statistics.Types (Estimator, Sample)++-- | A point and interval estimate computed via an 'Estimator'.+data Estimate = Estimate {+ estPoint :: {-# UNPACK #-} !Double+ -- ^ Point estimate.+ , estLowerBound :: {-# UNPACK #-} !Double+ -- ^ Lower bound of the estimate interval (i.e. the lower bound of+ -- the confidence interval).+ , estUpperBound :: {-# UNPACK #-} !Double+ -- ^ Upper bound of the estimate interval (i.e. the upper bound of+ -- the confidence interval).+ , estConfidenceLevel :: {-# UNPACK #-} !Double+ -- ^ Confidence level of the confidence intervals.+ } deriving (Eq, Show)++estimate :: Double -> Double -> Double -> Double -> Estimate+estimate pt lb ub cl =+ assert (lb <= ub) .+ assert (cl > 0 && cl < 1) $+ Estimate { estPoint = pt+ , estLowerBound = lb+ , estUpperBound = ub+ , estConfidenceLevel = cl+ }++data T = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double+infixl 2 :<++-- | Bias-corrected accelerated (BCA) bootstrap. This adjusts for both+-- bias and skewness in the resampled distribution.+bootstrapBCA :: Double -- ^ Confidence level+ -> Sample -- ^ Sample data+ -> [Estimator] -- ^ Estimators+ -> [Resample] -- ^ Resampled data+ -> [Estimate]+bootstrapBCA confidenceLevel sample =+ assert (confidenceLevel > 0 && confidenceLevel < 1)+ zipWith e+ where+ e est (Resample resample)+ | lengthU sample == 1 = estimate pt pt pt confidenceLevel+ | otherwise = + estimate pt (indexU resample lo) (indexU resample hi) confidenceLevel+ where+ pt = est sample+ lo = max (cumn a1) 0+ where a1 = bias + b1 / (1 - accel * b1)+ b1 = bias + z1+ hi = min (cumn a2) (ni - 1)+ where a2 = bias + b2 / (1 - accel * b2)+ b2 = bias - z1+ z1 = inverse standard ((1 - confidenceLevel) / 2)+ cumn = round . (*n) . cumulative standard+ bias = inverse standard (probN / n)+ where probN = fromIntegral . lengthU . filterU (<pt) $ resample+ ni = lengthU resample+ n = fromIntegral ni+ accel = sumCubes / (6 * (sumSquares ** 1.5))+ where (sumSquares :< sumCubes) = foldlU f (0 :< 0) jack+ f (s :< c) j = s + d2 :< c + d2 * d+ where d = jackMean - j+ d2 = d * d+ jackMean = mean jack+ jack = jackknife est sample++-- $references+--+-- * Davison, A.C; Hinkley, D.V. (1997) Bootstrap methods and their+-- application. <http://statwww.epfl.ch/davison/BMA/>
Statistics/Sample.hs view
@@ -12,8 +12,10 @@ module Statistics.Sample (+ -- * Types+ Sample -- * Statistics of location- mean+ , mean , harmonicMean , geometricMean @@ -42,28 +44,29 @@ -- | Arithmetic mean. This uses Welford's algorithm to provide -- numerical stability, using a single pass over the sample data. mean :: Sample -> Double-mean = fstT . foldlU k (T 0 0)- where- k (T m n) x = T m' n'- where m' = m + (x - m) / fromIntegral n'- n' = n + 1+mean = fini . foldlU go (T 0 0)+ where+ fini (T a _) = a+ go (T m n) x = T m' n'+ where m' = m + (x - m) / fromIntegral n'+ n' = n + 1 {-# INLINE mean #-} -- | Harmonic mean. This algorithm performs a single pass over the -- sample. harmonicMean :: Sample -> Double-harmonicMean xs = fromIntegral a / b+harmonicMean = fini . foldlU go (T 0 0) where- T b a = foldlU k (T 0 0) xs- k (T b a) n = T (b + (1/n)) (a+1)+ fini (T b a) = fromIntegral a / b+ go (T x y) n = T (x + (1/n)) (y+1) {-# INLINE harmonicMean #-} -- | Geometric mean of a sample containing no negative values. geometricMean :: Sample -> Double-geometricMean xs = p ** (1 / fromIntegral n)+geometricMean = fini . foldlU go (T 1 0) where- T p n = foldlU k (T 1 0) xs- k (T p n) a = T (p * a) (n + 1)+ fini (T p n) = p ** (1 / fromIntegral n)+ go (T p n) a = T (p * a) (n + 1) {-# INLINE geometricMean #-} -- $variance@@ -81,7 +84,7 @@ -- subject to stream fusion. robustVar :: Sample -> T-robustVar s = fini . foldlU go (T1 0 0 0) $ s+robustVar samp = fini . foldlU go (T1 0 0 0) $ samp where go (T1 n s c) x = T1 n' s' c' where n' = n + 1@@ -89,7 +92,7 @@ c' = c + d d = x - m fini (T1 n s c) = T (s - c ** (2 / fromIntegral n)) n- m = mean s+ m = mean samp -- | Maximum likelihood estimate of a sample's variance. variance :: Sample -> Double@@ -108,7 +111,7 @@ {-# INLINE varianceUnbiased #-} -- | Standard deviation. This is simply the square root of the--- maximum likelihood estimate of the variance. +-- maximum likelihood estimate of the variance. stdDev :: Sample -> Double stdDev = sqrt . varianceUnbiased @@ -149,7 +152,7 @@ {-# INLINE fastVarianceUnbiased #-} -- | Standard deviation. This is simply the square root of the--- maximum likelihood estimate of the variance. +-- maximum likelihood estimate of the variance. fastStdDev :: Sample -> Double fastStdDev = sqrt . fastVariance {-# INLINE fastStdDev #-}@@ -161,9 +164,6 @@ data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int data T1 = T1 {-# UNPACK #-}!Int {-# UNPACK #-}!Double {-# UNPACK #-}!Double--fstT :: T -> Double-fstT (T a _) = a {-
Statistics/Types.hs view
@@ -12,10 +12,18 @@ module Statistics.Types ( Sample+ , Estimator , Weights ) where import Data.Array.Vector (UArr) +-- | Sample data. type Sample = UArr Double++-- | A function that estimates a property of a sample, such as its+-- 'mean'.+type Estimator = Sample -> Double++-- | Weights for affecting the importance of elements of a sample. type Weights = UArr Double
statistics.cabal view
@@ -1,7 +1,20 @@ name: statistics-version: 0.1-synopsis: A library of statistical types, data, and functions.-description: A library of statistical types, data, and functions.+version: 0.2+synopsis: A library of statistical types, data, and functions+description:+ This library provides a number of common functions and types useful+ in statistics. Our focus is on high performance, numerical+ 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:+ .+ Working with widely used discrete and continuous probability+ distributions. (There are dozens of exotic distributions in use; we+ focus on the most common.)+ .+ Computing with sample data: quantile estimation, kernel density+ estimation, bootstrap methods, and autocorrelation analysis. license: BSD3 license-file: LICENSE homepage: http://darcs.serpentine.com/statistics@@ -15,20 +28,30 @@ library exposed-modules:+ Statistics.Autocorrelation+ Statistics.Constants Statistics.Distribution+ Statistics.Distribution.Binomial+ Statistics.Distribution.Gamma+ Statistics.Distribution.Geometric+ Statistics.Distribution.Exponential+ Statistics.Distribution.Hypergeometric Statistics.Distribution.Normal+ Statistics.Distribution.Poisson Statistics.Function Statistics.KernelDensity+ Statistics.Math Statistics.Quantile+ Statistics.Resampling+ Statistics.Resampling.Bootstrap Statistics.Sample Statistics.Types- other-modules:- Statistics.Constants build-depends: base < 5, erf,+ mersenne-random, uvector >= 0.1.0.4,- uvector-algorithms+ uvector-algorithms >= 0.2 if impl(ghc >= 6.10) build-depends: base >= 4