packages feed

statistics (empty) → 0.1

raw patch · 12 files changed

+846/−0 lines, 12 filesdep +basedep +erfdep +uvectorsetup-changed

Dependencies added: base, erf, uvector, uvector-algorithms

Files

+ LICENSE view
@@ -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.
+ README view
@@ -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 <bos@serpentine.com>
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ Statistics/Constants.hs view
@@ -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 #-}
+ Statistics/Distribution.hs view
@@ -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<X)@.+    cumulative  :: d -> Double -> Double++    -- | Inverse of the cumulative distribution function.  The value+    -- @X@ for which @P(x<X)@.+    inverse     :: d -> 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
+ Statistics/Distribution/Normal.hs view
@@ -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
+ Statistics/Function.hs view
@@ -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 #-}
+ Statistics/KernelDensity.hs view
@@ -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")
+ Statistics/Quantile.hs view
@@ -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/.+--   <http://mathworld.wolfram.com/Quantile.html>+--+-- * Hyndman, R.J.; Fan, Y. (1996) Sample quantiles in statistical+--   packages. /American Statistician/+--   50(4):361&#8211;365. <http://www.jstor.org/stable/2684934>+
+ Statistics/Sample.hs view
@@ -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&#8212;and hence the standard deviation&#8212;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. <ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>+--+-- * 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&#8211;420. <http://www.jstor.org/stable/1266577>+--+-- * West, D.H.D. (1979) Updating mean and variance estimates: an+--   improved method. /Communications of the ACM/+--   22(9):532&#8211;535. <http://doi.acm.org/10.1145/359146.359153>
+ Statistics/Types.hs view
@@ -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
+ statistics.cabal view
@@ -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 <bos@serpentine.com>+maintainer:     Bryan O'Sullivan <bos@serpentine.com>+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