diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,26 @@ +Copyright (c) 2009, Bryan O'Sullivan +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,20 @@ +Statistics: efficient, general purpose statistics +------------------------------------------------- + +This package provides the Statistics module, a Haskell library for +working with statistical data in a space- and time-efficient way. + +Where possible, we give citations and computational complexity +estimates for the algorithms used. + + +Source code +----------- + +darcs get http://darcs.serpentine.com/statistics + + +Authors +------- + +Bryan O'Sullivan diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/Statistics/Constants.hs b/Statistics/Constants.hs new file mode 100644 --- /dev/null +++ b/Statistics/Constants.hs @@ -0,0 +1,44 @@ +-- | +-- Module : Statistics.Constants +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Constant values common to much statistics code. + +module Statistics.Constants + ( + m_huge + , m_1_sqrt_2 + , m_2_sqrt_pi + , m_sqrt_2 + , m_sqrt_2_pi + ) where + +-- | A very large number. +m_huge :: Double +m_huge = 1.797693e308 +{-# INLINE m_huge #-} + +-- | @sqrt 2@ +m_sqrt_2 :: Double +m_sqrt_2 = 1.414213562373095145474621858739 +{-# INLINE m_sqrt_2 #-} + +-- | @sqrt (2 * pi)@ +m_sqrt_2_pi :: Double +m_sqrt_2_pi = 2.506628274631000241612355239340 +{-# INLINE m_sqrt_2_pi #-} + +-- | @2 / sqrt pi@ +m_2_sqrt_pi :: Double +m_2_sqrt_pi = 1.128379167095512558560699289956 +{-# INLINE m_2_sqrt_pi #-} + +-- | @1 / sqrt 2@ +m_1_sqrt_2 :: Double +m_1_sqrt_2 = 0.707106781186547461715008466854 +{-# INLINE m_1_sqrt_2 #-} diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution.hs @@ -0,0 +1,61 @@ +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} +-- | +-- Module : Statistics.Distribution +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Types and functions common to many probability distributions. + +module Statistics.Distribution + ( + Distribution(..) + , 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)@. + 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@. +-- +-- 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@. +findRoot :: Distribution d => d + -> Double -- ^ Probability @p@ + -> Double -- ^ Initial guess + -> Double -- ^ Lower bound on interval + -> Double -- ^ Upper bound on interval + -> Double +findRoot d prob = loop 0 1 + where + loop !(i::Int) !dx !x !lo !hi + | abs dx <= accuracy || i >= maxIters = x + | otherwise = loop (i+1) dx'' x'' lo' hi' + where + err = cumulative d x - prob + (lo',hi') | err < 0 = (x, hi) + | otherwise = (lo, x) + pdf = probability d x + (dx',x') | pdf /= 0 = (err / pdf, x - dx) + | otherwise = (dx, x) + (dx'',x'') + | x' < lo' || x' > hi' || pdf == 0 = (x'-x, (lo + hi) / 2) + | otherwise = (dx', x') + accuracy = 1e-15 + maxIters = 150 diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs new file mode 100644 --- /dev/null +++ b/Statistics/Distribution/Normal.hs @@ -0,0 +1,73 @@ +-- | +-- Module : Statistics.Normal +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- The normal distribution. + +module Statistics.Distribution.Normal + ( + NormalDistribution + , fromParams + , fromSample + , standard + ) where + +import Control.Exception (assert) +import Data.Number.Erf (erfc) +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 { + mean :: {-# UNPACK #-} !Double + , variance :: {-# UNPACK #-} !Double + , pdfDenom :: {-# UNPACK #-} !Double + , cdfDenom :: {-# UNPACK #-} !Double + } deriving (Eq, Ord, Read, Show) + +instance D.Distribution NormalDistribution where + probability = probability + cumulative = cumulative + inverse = inverse + +standard :: NormalDistribution +standard = NormalDistribution { + mean = 0.0 + , variance = 1.0 + , cdfDenom = m_sqrt_2 + , pdfDenom = m_sqrt_2_pi + } + +fromParams :: Double -> Double -> NormalDistribution +fromParams m v = assert (v > 0) $ + NormalDistribution { + mean = m + , variance = v + , cdfDenom = m_sqrt_2 * sv + , pdfDenom = m_sqrt_2_pi * 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 + where xm = x - mean d + +cumulative :: NormalDistribution -> Double -> Double +cumulative d x = erfc (-(x-mean d) / cdfDenom d) / 2 + +inverse :: NormalDistribution -> Double -> Double +inverse d p + | p == 0 = -m_huge + | p == 1 = m_huge + | p == 0.5 = mean d + | otherwise = x * sqrt (variance d) + mean d + where x = D.findRoot standard p 0 (-100) 100 diff --git a/Statistics/Function.hs b/Statistics/Function.hs new file mode 100644 --- /dev/null +++ b/Statistics/Function.hs @@ -0,0 +1,46 @@ +{-# LANGUAGE TypeOperators #-} +-- | +-- Module : Statistics.Quantile +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Functions for computing quantiles. + +module Statistics.Function + ( + minMax + , sort + , partialSort + ) where + +import Data.Array.Vector.Algorithms.Immutable (apply) +import Data.Array.Vector ((:*:)(..), UA, UArr, foldlU) +import qualified Data.Array.Vector.Algorithms.Intro as I + +-- | Sort. +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 +-- at the front. +partialSort :: (UA e, Ord e) => + Int -- ^ The number @k@ of least elements + -> UArr e + -> UArr e +partialSort k = apply (\a -> I.partialSort a k) +{-# INLINE partialSort #-} + +data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double + +-- | Compute the minimum and maximum of an array in one pass. +minMax :: UArr Double -> Double :*: Double +minMax = fini . foldlU go (MM (1/0) (-1/0)) + where + go (MM lo hi) k = MM (min lo k) (max hi k) + fini (MM lo hi) = lo :*: hi +{-# INLINE minMax #-} diff --git a/Statistics/KernelDensity.hs b/Statistics/KernelDensity.hs new file mode 100644 --- /dev/null +++ b/Statistics/KernelDensity.hs @@ -0,0 +1,161 @@ +-- | +-- Module : Statistics.KernelDensity +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Kernel density estimation code, providing non-parametric ways to +-- estimate the probability density function of a sample. + +module Statistics.KernelDensity + ( + -- * Simple entry points + epanechnikovPDF + , gaussianPDF + -- * Building blocks + -- These functions may be useful if you need to construct a kernel + -- density function estimator other than the ones provided in this + -- module. + + -- ** Choosing points from a sample + , Points(..) + , choosePoints + -- ** Bandwidth estimation + , Bandwidth + , bandwidth + , epanechnikovBW + , gaussianBW + -- ** Kernels + , Kernel + , epanechnikovKernel + , gaussianKernel + -- ** Low-level estimation + , estimatePDF + , simplePDF + ) where + +import Data.Array.Vector ((:*:)(..), UArr, enumFromToU, lengthU, mapU, sumU) +import Statistics.Function (minMax) +import Statistics.Sample (stdDev) +import Statistics.Constants (m_1_sqrt_2, m_2_sqrt_pi) +import Statistics.Types (Sample) + +-- | Points from the range of a 'Sample'. +newtype Points = Points { + fromPoints :: UArr Double + } deriving (Eq, Show) + +-- | Bandwidth estimator for an Epanechnikov kernel. +epanechnikovBW :: Double -> Bandwidth +epanechnikovBW n = (80 / (n * m_2_sqrt_pi)) ** 0.2 + +-- | Bandwidth estimator for a Gaussian kernel. +gaussianBW :: Double -> Bandwidth +gaussianBW n = (4 / (n * 3)) ** 0.2 + +-- | The width of the convolution kernel used. +type Bandwidth = Double + +-- | Compute the optimal bandwidth from the observed data for the given +-- kernel. +bandwidth :: (Double -> Bandwidth) + -> Sample + -> Bandwidth +bandwidth kern values = stdDev values * kern (fromIntegral $ lengthU values) + +-- | Choose a uniform range of points at which to estimate a sample's +-- probability density function. +-- +-- If you are using a Gaussian kernel, multiply the sample's bandwidth +-- by 3 before passing it to this function. +-- +-- If this function is passed an empty vector, it returns values of +-- positive and negative infinity. +choosePoints :: Int -- ^ Number of points to select, /n/ + -> Double -- ^ Sample bandwidth, /h/ + -> Sample -- ^ Input data + -> Points +choosePoints n h sample = Points . mapU f $ enumFromToU 0 n' + where lo = a - h + hi = z + h + a :*: z = minMax sample + d = (hi - lo) / fromIntegral n' + f i = lo + fromIntegral i * d + 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 + -> Double + -> Double + -> Double + +-- | Epanechnikov kernel for probability density function estimation. +epanechnikovKernel :: Kernel +epanechnikovKernel f h p v + | abs u <= 1 = f * (1 - u * u) + | otherwise = 0 + where u = (v - p) / (h * 0.75) + +-- | Gaussian kernel for probability density function estimation. +gaussianKernel :: Kernel +gaussianKernel f h p v = exp (-0.5 * u * u) * g + where u = (v - p) / h + g = f * m_2_sqrt_pi * m_1_sqrt_2 + +-- | Kernel density estimator, providing a non-parametric way of +-- estimating the PDF of a random variable. +estimatePDF :: Kernel -- ^ Kernel function + -> Bandwidth -- ^ Bandwidth, /h/ + -> Sample -- ^ Sample data + -> Points -- ^ Points at which to estimate + -> UArr Double +estimatePDF kernel h sample + | n < 2 = errorShort "estimatePDF" + | otherwise = mapU k . fromPoints + where + k p = sumU . mapU (kernel f h p) $ sample + f = 1 / (h * fromIntegral n) + n = lengthU sample +{-# INLINE estimatePDF #-} + +-- | A helper for creating a simple kernel density estimation function +-- with automatically chosen bandwidth and estimation points. +simplePDF :: (Double -> Double) -- ^ Bandwidth function + -> Kernel -- ^ Kernel function + -> Double -- ^ Bandwidth scaling factor (3 for a Gaussian kernel, 1 for all others) + -> Int -- ^ Number of points at which to estimate + -> Sample -- ^ Sample data + -> (Points, UArr Double) +simplePDF fbw fpdf k numPoints sample = + (points, estimatePDF fpdf bw sample points) + where points = choosePoints numPoints (bw*k) sample + bw = bandwidth fbw sample +{-# INLINE simplePDF #-} + +-- | Simple Epanechnikov kernel density estimator. Returns the +-- uniformly spaced points from the sample range at which the density +-- function was estimated, and the estimates at those points. +epanechnikovPDF :: Int -- ^ Number of points at which to estimate + -> Sample + -> (Points, UArr Double) +epanechnikovPDF = simplePDF epanechnikovBW epanechnikovKernel 1 + +-- | Simple Gaussian kernel density estimator. Returns the uniformly +-- spaced points from the sample range at which the density function +-- was estimated, and the estimates at those points. +gaussianPDF :: Int -- ^ Number of points at which to estimate + -> Sample + -> (Points, UArr Double) +gaussianPDF = simplePDF gaussianBW gaussianKernel 3 + +errorShort :: String -> a +errorShort func = error ("Statistics.KernelDensity." ++ func ++ + ": at least two points required") diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs new file mode 100644 --- /dev/null +++ b/Statistics/Quantile.hs @@ -0,0 +1,143 @@ +{-# LANGUAGE TypeOperators #-} +-- | +-- Module : Statistics.Quantile +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Functions for approximating quantiles. + +module Statistics.Quantile + ( + -- * Types + ContParam(..) + + -- * Quantile estimation functions + , weightedAvg + , continuousBy + + -- * Parameters for the continuous sample method + , cadpw + , hazen + , s + , spss + , medianUnbiased + , normalUnbiased + + -- * References + -- $references + ) where + +import Control.Exception (assert) +import Data.Array.Vector (allU, indexU, lengthU) +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 + -> Double +weightedAvg k q x = + assert (q >= 2) . + assert (k >= 0) . + assert (k < q) . + assert (allU (not . isNaN) x) $ + xj + g * (xj1 - xj) + where + j = floor idx + idx = fromIntegral (lengthU x - 1) * fromIntegral k / fromIntegral q + g = idx - fromIntegral j + xj = indexU sx j + xj1 = indexU sx (j+1) + sx = partialSort (j+2) x +{-# INLINE weightedAvg #-} + +-- | Parameters @a@ and @b@ to the 'quantileBy' 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 + -> Double +continuousBy (ContParam a b) k q x = + assert (q >= 2) . + assert (k >= 0) . + assert (k <= q) . + assert (allU (not . isNaN) x) $ + (1-h) * item (j-1) + h * item j + where + j = floor (t + eps) + t = a + p * (fromIntegral n + 1 - a - b) + p = fromIntegral k / fromIntegral q + h | abs r < eps = 0 + | otherwise = r + where r = t - fromIntegral j + eps = 8.881784e-16 + n = lengthU x + item m = indexU sx $ bracket m + 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. +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 +-- 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. +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. +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 +-- Mathematica. +medianUnbiased :: ContParam +medianUnbiased = ContParam third third + where third = 1/3 +{-# INLINE medianUnbiased #-} + +-- | 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. +normalUnbiased :: ContParam +normalUnbiased = ContParam ta ta + where ta = 3/8 +{-# INLINE normalUnbiased #-} + +-- $references +-- +-- * Weisstein, E.W. Quantile. /MathWorld/. +-- +-- +-- * Hyndman, R.J.; Fan, Y. (1996) Sample quantiles in statistical +-- packages. /American Statistician/ +-- 50(4):361–365. + diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs new file mode 100644 --- /dev/null +++ b/Statistics/Sample.hs @@ -0,0 +1,207 @@ +-- | +-- Module : Statistics.Sample +-- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Commonly used sample statistics, also known as descriptive +-- statistics. + +module Statistics.Sample + ( + -- * Statistics of location + mean + , harmonicMean + , geometricMean + + -- * Statistics of dispersion + -- $variance + + -- ** Two-pass functions (numerically robust) + -- $robust + , variance + , varianceUnbiased + , stdDev + + -- ** Single-pass functions (faster, less safe) + -- $cancellation + , fastVariance + , fastVarianceUnbiased + , fastStdDev + + -- * References + -- $references + ) where + +import Data.Array.Vector (foldlU) +import Statistics.Types (Sample) + +-- | 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 +{-# INLINE mean #-} + +-- | Harmonic mean. This algorithm performs a single pass over the +-- sample. +harmonicMean :: Sample -> Double +harmonicMean xs = fromIntegral a / b + where + T b a = foldlU k (T 0 0) xs + k (T b a) n = T (b + (1/n)) (a+1) +{-# INLINE harmonicMean #-} + +-- | Geometric mean of a sample containing no negative values. +geometricMean :: Sample -> Double +geometricMean xs = p ** (1 / fromIntegral n) + where + T p n = foldlU k (T 1 0) xs + k (T p n) a = T (p * a) (n + 1) +{-# INLINE geometricMean #-} + +-- $variance +-- +-- The variance—and hence the standard deviation—of a +-- sample of fewer than two elements are both defined to be zero. + +-- $robust +-- +-- These functions use the compensated summation algorithm of Chan et +-- al. for numerical robustness, but require two passes over the +-- sample data as a result. +-- +-- Because of the need for two passes, these functions are /not/ +-- subject to stream fusion. + +robustVar :: Sample -> T +robustVar s = fini . foldlU go (T1 0 0 0) $ s + where + go (T1 n s c) x = T1 n' s' c' + where n' = n + 1 + s' = s + d * d + c' = c + d + d = x - m + fini (T1 n s c) = T (s - c ** (2 / fromIntegral n)) n + m = mean s + +-- | Maximum likelihood estimate of a sample's variance. +variance :: Sample -> Double +variance = fini . robustVar + where fini (T v n) + | n > 1 = v / fromIntegral n + | otherwise = 0 +{-# INLINE variance #-} + +-- | Unbiased estimate of a sample's variance. +varianceUnbiased :: Sample -> Double +varianceUnbiased = fini . robustVar + where fini (T v n) + | n > 1 = v / fromIntegral (n-1) + | otherwise = 0 +{-# INLINE varianceUnbiased #-} + +-- | Standard deviation. This is simply the square root of the +-- maximum likelihood estimate of the variance. +stdDev :: Sample -> Double +stdDev = sqrt . varianceUnbiased + +-- $cancellation +-- +-- The functions prefixed with the name @fast@ below perform a single +-- pass over the sample data using Knuth's algorithm. They usually +-- work well, but see below for caveats. These functions are subject +-- to array fusion. +-- +-- /Note/: in cases where most sample data is close to the sample's +-- mean, Knuth's algorithm gives inaccurate results due to +-- catastrophic cancellation. + +fastVar :: Sample -> T1 +fastVar = foldlU go (T1 0 0 0) + where + go (T1 n m s) x = T1 n' m' s' + where n' = n + 1 + m' = m + d / fromIntegral n' + s' = s + d * (x - m') + d = x - m + +-- | Maximum likelihood estimate of a sample's variance. +fastVariance :: Sample -> Double +fastVariance = fini . fastVar + where fini (T1 n _m s) + | n > 1 = s / fromIntegral n + | otherwise = 0 +{-# INLINE fastVariance #-} + +-- | Unbiased estimate of a sample's variance. +fastVarianceUnbiased :: Sample -> Double +fastVarianceUnbiased = fini . fastVar + where fini (T1 n _m s) + | n > 1 = s / fromIntegral (n - 1) + | otherwise = 0 +{-# INLINE fastVarianceUnbiased #-} + +-- | Standard deviation. This is simply the square root of the +-- maximum likelihood estimate of the variance. +fastStdDev :: Sample -> Double +fastStdDev = sqrt . fastVariance +{-# INLINE fastStdDev #-} + +------------------------------------------------------------------------ +-- Helper code. Monomorphic unpacked accumulators. + +-- don't support polymorphism, as we can't get unboxed returns if we use it. +data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int + +data T1 = T1 {-# UNPACK #-}!Int {-# UNPACK #-}!Double {-# UNPACK #-}!Double + +fstT :: T -> Double +fstT (T a _) = a + +{- + +Consider this core: + +with data T a = T !a !Int + +$wfold :: Double# + -> Int# + -> Int# + -> (# Double, Int# #) + +and without, + +$wfold :: Double# + -> Int# + -> Int# + -> (# Double#, Int# #) + +yielding to boxed returns and heap checks. + +-} + +-- $references +-- +-- * Chan, T. F.; Golub, G.H.; LeVeque, R.J. (1979) Updating formulae +-- and a pairwise algorithm for computing sample +-- variances. Technical Report STAN-CS-79-773, Department of +-- Computer Science, Stanford +-- University. +-- +-- * Knuth, D.E. (1998) The art of computer programming, volume 2: +-- seminumerical algorithms, 3rd ed., p. 232. +-- +-- * Welford, B.P. (1962) Note on a method for calculating corrected +-- sums of squares and products. /Technometrics/ +-- 4(3):419–420. +-- +-- * West, D.H.D. (1979) Updating mean and variance estimates: an +-- improved method. /Communications of the ACM/ +-- 22(9):532–535. diff --git a/Statistics/Types.hs b/Statistics/Types.hs new file mode 100644 --- /dev/null +++ b/Statistics/Types.hs @@ -0,0 +1,21 @@ +-- | +-- Module : Statistics.Types +-- Copyright : (c) 2009 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Types for working with statistics. + +module Statistics.Types + ( + Sample + , Weights + ) where + +import Data.Array.Vector (UArr) + +type Sample = UArr Double +type Weights = UArr Double diff --git a/statistics.cabal b/statistics.cabal new file mode 100644 --- /dev/null +++ b/statistics.cabal @@ -0,0 +1,41 @@ +name: statistics +version: 0.1 +synopsis: A library of statistical types, data, and functions. +description: A library of statistical types, data, and functions. +license: BSD3 +license-file: LICENSE +homepage: http://darcs.serpentine.com/statistics +author: Bryan O'Sullivan +maintainer: Bryan O'Sullivan +copyright: 2009 Bryan O'Sullivan +category: Math, Statistics +build-type: Simple +cabal-version: >= 1.2 +extra-source-files: README + +library + exposed-modules: + Statistics.Distribution + Statistics.Distribution.Normal + Statistics.Function + Statistics.KernelDensity + Statistics.Quantile + Statistics.Sample + Statistics.Types + other-modules: + Statistics.Constants + build-depends: + base < 5, + erf, + uvector >= 0.1.0.4, + uvector-algorithms + if impl(ghc >= 6.10) + build-depends: + base >= 4 + + -- gather extensive profiling data for now + ghc-prof-options: -auto-all + + ghc-options: -Wall -funbox-strict-fields -O2 + if impl(ghc >= 6.8) + ghc-options: -fwarn-tabs