diff --git a/Statistics/Autocorrelation.hs b/Statistics/Autocorrelation.hs new file mode 100644 --- /dev/null +++ b/Statistics/Autocorrelation.hs @@ -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 diff --git a/Statistics/Constants.hs b/Statistics/Constants.hs --- a/Statistics/Constants.hs +++ b/Statistics/Constants.hs @@ -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) diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs --- a/Statistics/Distribution.hs +++ b/Statistics/Distribution.hs @@ -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 Double -> Double -- | Inverse of the cumulative distribution function. The value - -- @X@ for which @P(x 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 diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Binomial.hs @@ -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 #-} diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Exponential.hs @@ -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 #-} diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Gamma.hs @@ -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 #-} diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Geometric.hs @@ -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 #-} diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Hypergeometric.hs @@ -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 #-} diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs --- a/Statistics/Distribution/Normal.hs +++ b/Statistics/Distribution/Normal.hs @@ -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 diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Poisson.hs @@ -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 #-} diff --git a/Statistics/Function.hs b/Statistics/Function.hs --- a/Statistics/Function.hs +++ b/Statistics/Function.hs @@ -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) diff --git a/Statistics/KernelDensity.hs b/Statistics/KernelDensity.hs --- a/Statistics/KernelDensity.hs +++ b/Statistics/KernelDensity.hs @@ -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 diff --git a/Statistics/Math.hs b/Statistics/Math.hs new file mode 100644 --- /dev/null +++ b/Statistics/Math.hs @@ -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. +-- +-- * 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. +-- +-- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete +-- gamma integral. /Applied Statistics/ +-- 37(3):466–473. diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs --- a/Statistics/Quantile.hs +++ b/Statistics/Quantile.hs @@ -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. - diff --git a/Statistics/Resampling.hs b/Statistics/Resampling.hs new file mode 100644 --- /dev/null +++ b/Statistics/Resampling.hs @@ -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 diff --git a/Statistics/Resampling/Bootstrap.hs b/Statistics/Resampling/Bootstrap.hs new file mode 100644 --- /dev/null +++ b/Statistics/Resampling/Bootstrap.hs @@ -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 ( diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs --- a/Statistics/Sample.hs +++ b/Statistics/Sample.hs @@ -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 {- diff --git a/Statistics/Types.hs b/Statistics/Types.hs --- a/Statistics/Types.hs +++ b/Statistics/Types.hs @@ -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 diff --git a/statistics.cabal b/statistics.cabal --- a/statistics.cabal +++ b/statistics.cabal @@ -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