statistics 0.2.1 → 0.2.2

raw patch · 7 files changed

+405/−20 lines, 7 files

Files

Statistics/Distribution/Exponential.hs +8/−6
Statistics/Distribution/Geometric.hs +8/−4
Statistics/Internal.hs +41/−0
Statistics/Quantile.hs +39/−6
Statistics/Sample.hs +93/−3
Statistics/Sample/Powers.hs +212/−0
statistics.cabal +4/−1

Statistics/Distribution/Exponential.hs view

@@ -8,10 +8,10 @@ -- 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/.+-- The exponential distribution.  This is the continunous probability+-- distribution of the times between events in a poisson process, in+-- which events occur continuously and independently at a constant+-- average rate.  module Statistics.Distribution.Exponential     (@@ -19,6 +19,8 @@     -- * Constructors     , fromLambda     , fromSample+    -- * Accessors+    , edLambda     ) where  import Data.Typeable (Typeable)@@ -39,11 +41,11 @@     {-# INLINE inverse #-}  instance D.Variance ExponentialDistribution where-    variance (ED l) = l * l+    variance (ED l) = 1 / (l * l)     {-# INLINE variance #-}  instance D.Mean ExponentialDistribution where-    mean = edLambda+    mean (ED l) = 1 / l     {-# INLINE mean #-}  fromLambda :: Double            -- ^ &#955; (scale) parameter.

Statistics/Distribution/Geometric.hs view

@@ -8,10 +8,14 @@ -- 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.+-- The Geometric distribution. This is the probability distribution of+-- the number of Bernoulli trials needed to get one success, supported+-- on the set [1,2..].+--+-- This distribution is sometimes referred to as the /shifted/+-- geometric distribution, to distinguish it from a variant measuring+-- the number of failures before the first success, defined over the+-- set [0,1..].  module Statistics.Distribution.Geometric     (

+ Statistics/Internal.hs view

@@ -0,0 +1,41 @@+{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-}+-- |+-- Module    : Statistics.Internal+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Scary internal functions.++module Statistics.Internal+    (+      inlinePerformIO+    ) where++#if __GLASGOW_HASKELL__ >= 611+import GHC.IO (IO(IO))+#else+import GHC.IOBase (IO(IO))+#endif+import GHC.Base (realWorld#)+#if !defined(__GLASGOW_HASKELL__)+import System.IO.Unsafe (unsafePerformIO)+#endif++-- Lifted from Data.ByteString.Internal so we don't introduce an+-- otherwise unnecessary dependency on the bytestring package.++-- | Just like unsafePerformIO, but we inline it. Big performance+-- gains as it exposes lots of things to further inlining. /Very+-- unsafe/. In particular, you should do no memory allocation inside+-- an 'inlinePerformIO' block. On Hugs this is just @unsafePerformIO@.+{-# INLINE inlinePerformIO #-}+inlinePerformIO :: IO a -> a+#if defined(__GLASGOW_HASKELL__)+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r+#else+inlinePerformIO = unsafePerformIO+#endif

Statistics/Quantile.hs view

@@ -23,6 +23,7 @@       weightedAvg     , ContParam(..)     , continuousBy+    , midspread      -- * Parameters for the continuous sample method     , cadpw@@ -42,8 +43,8 @@ import Statistics.Function (partialSort) import Statistics.Types (Sample) --- | Estimate the /k/th /q/-quantile of a sample, using the weighted--- average method.+-- | O(/n/ log /n/). 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.@@ -66,10 +67,10 @@ -- | Parameters /a/ and /b/ to the 'continuousBy' function. data ContParam = ContParam {-# UNPACK #-} !Double {-# UNPACK #-} !Double --- | 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.+-- | O(/n/ log /n/). 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.@@ -94,6 +95,38 @@     sx              = partialSort (bracket j + 1) x     bracket m       = min (max m 0) (n - 1) {-# INLINE continuousBy #-}++-- | O(/n/ log /n/). Estimate the range between /q/-quantiles 1 and+-- /q/-1 of a sample /x/, using the continuous sample method with the+-- given parameters.+--+-- For instance, the interquartile range (IQR) can be estimated as+-- follows:+--+-- > midspread medianUnbiased 4 (toU [1,1,2,2,3])+-- > ==> 1.333333+midspread :: ContParam       -- ^ Parameters /a/ and /b/.+          -> Int             -- ^ /q/, the number of quantiles.+          -> Sample          -- ^ /x/, the sample data.+          -> Double+midspread (ContParam a b) k x =+    assert (allU (not . isNaN) x) .+    assert (k > 0) $+    quantile (1-frac) - quantile frac+  where+    quantile i        = (1-h i) * item (j i-1) + h i * item (j i)+    j i               = floor (t i + eps) :: Int+    t i               = a + i * (fromIntegral n + 1 - a - b)+    h i | abs r < eps = 0+        | otherwise   = r+        where r       = t i - fromIntegral (j i)+    eps               = m_epsilon * 4+    n                 = lengthU x+    item              = indexU sx . bracket+    sx                = partialSort (bracket (j (1-frac)) + 1) x+    bracket m         = min (max m 0) (n - 1)+    frac              = 1 / fromIntegral k+{-# INLINE midspread #-}  -- | California Department of Public Works definition, /a/=0, /b/=1. -- Gives a linear interpolation of the empirical CDF.  This

Statistics/Sample.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE TypeOperators #-} -- | -- Module    : Statistics.Sample -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan@@ -14,6 +15,9 @@     (     -- * Types       Sample+    -- * Descriptive functions+    , range+     -- * Statistics of location     , mean     , harmonicMean@@ -22,6 +26,12 @@     -- * Statistics of dispersion     -- $variance +    -- ** Functions over central moments+    , centralMoment+    , centralMoments+    , skewness+    , kurtosis+     -- ** Two-pass functions (numerically robust)     -- $robust     , variance@@ -38,9 +48,15 @@     -- $references     ) where -import Data.Array.Vector (foldlU, lengthU)+import Data.Array.Vector+import Statistics.Function (minMax) import Statistics.Types (Sample) +range :: Sample -> Double+range s = hi - lo+    where hi :*: lo = minMax s+{-# INLINE range #-}+ -- | Arithmetic mean.  This uses Welford's algorithm to provide -- numerical stability, using a single pass over the sample data. mean :: Sample -> Double@@ -69,6 +85,78 @@     go (T p n) a = T (p * a) (n + 1) {-# INLINE geometricMean #-} +-- | Compute the /k/th central moment of a sample.+--+-- This function performs two passes over the sample, so is not subject+-- to stream fusion.+centralMoment :: Int -> Sample -> Double+centralMoment a xs+    | a < 0  = error "Statistics.Sample.centralMoment: negative input"+    | a == 0 = 1+    | a == 1 = 0+    | otherwise = sumU (mapU go xs) / fromIntegral (lengthU xs)+  where+    go x = (x-m) ^ a+    m    = mean xs+{-# INLINE centralMoment #-}++-- | Compute the /k/th and /j/th central moments of a sample.+--+-- This function performs two passes over the sample, so is not subject+-- to stream fusion.+--+-- For samples containing many values very close to the mean, this+-- function is subject to inaccuracy due to catastrophic cancellation.+centralMoments :: Int -> Int -> Sample -> Double :*: Double+centralMoments a b xs+    | a < 2 || b < 2 = centralMoment a xs :*: centralMoment b xs+    | otherwise      = fini . foldlU go (V 0 0) $ xs+  where go (V i j) x = V (i + d^a) (j + d^b)+            where d  = x - m+        fini (V i j) = i / n :*: j / n+        m            = mean xs+        n            = fromIntegral (lengthU xs)+{-# INLINE centralMoments #-}++-- | Compute the skewness of a sample. This is a measure of the+-- asymmetry of its distribution.+--+-- A sample with negative skew is said to be /left-skewed/.  Most of+-- its mass is on the right of the distribution, with the tail on the+-- left.+--+-- > skewness . powers 3 $ toU [1,100,101,102,103]+-- > ==> -1.497681449918257+--+-- A sample with positive skew is said to be /right-skewed/.+--+-- > skewness . powers 5 $ toU [1,2,3,4,100]+-- > ==> 1.4975367033335198+--+-- A sample's skewness is not defined if its 'variance' is zero.+--+-- This function performs two passes over the sample, so is not subject+-- to stream fusion.+skewness :: Sample -> Double+skewness xs = c3 * c2 ** (-1.5)+    where c3 :*: c2 = centralMoments 3 2 xs+{-# INLINE skewness #-}++-- | Compute the excess kurtosis of a sample.  This is a measure of+-- the \"peakedness\" of its distribution.  A high kurtosis indicates+-- that more of the sample's variance is due to infrequent severe+-- deviations, rather than more frequent modest deviations.+--+-- A sample's excess kurtosis is not defined if its 'variance' is+-- zero.+--+-- This function performs two passes over the sample, so is not subject+-- to stream fusion.+kurtosis :: Sample -> Double+kurtosis xs = c4 / (c2 * c2) - 3+    where c4 :*: c2 = centralMoments 4 2 xs+{-# INLINE kurtosis #-}+ -- $variance -- -- The variance&#8212;and hence the standard deviation&#8212;of a@@ -94,7 +182,8 @@     n            = lengthU samp     m            = mean samp --- | Maximum likelihood estimate of a sample's variance.+-- | Maximum likelihood estimate of a sample's variance.  Also known+-- as the population variance, where the denominator is /n/. variance :: Sample -> Double variance = fini . robustVar   where fini (T v n)@@ -102,7 +191,8 @@           | otherwise = 0 {-# INLINE variance #-} --- | Unbiased estimate of a sample's variance.+-- | Unbiased estimate of a sample's variance.  Also known as the+-- sample variance, where the denominator is /n/-1. varianceUnbiased :: Sample -> Double varianceUnbiased = fini . robustVar   where fini (T v n)

+ Statistics/Sample/Powers.hs view

@@ -0,0 +1,212 @@+{-# LANGUAGE BangPatterns, TypeOperators #-}+-- |+-- Module    : Statistics.Sample.Powers+-- Copyright : (c) 2009 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Very fast statistics over simple powers of a sample.  These can all+-- be computed efficiently in just a single pass over a sample, with+-- that pass subject to stream fusion.+--+-- The tradeoff is that some of these functions are less numerically+-- robust than their counterparts in the 'Statistics.Sample' module.+-- Where this is the case, the alternatives are noted.++module Statistics.Sample.Powers+    (+    -- * Types+      Sample+    , Powers++    -- * Constructor+    , powers++    -- * Descriptive functions+    , order+    , count+    , sum++    -- * Statistics of location+    , mean++    -- * Statistics of dispersion+    , variance+    , stdDev+    , varianceUnbiased++    -- * Functions over central moments+    , centralMoment+    , skewness+    , kurtosis++    -- * References+    -- $references+    ) where++import Control.Monad.ST (unsafeSTToIO)+import Data.Array.Vector+import Prelude hiding (sum)+import Statistics.Internal (inlinePerformIO)+import Statistics.Math (choose)+import Statistics.Types (Sample)+import System.IO.Unsafe (unsafePerformIO)++newtype Powers = Powers (UArr Double)+    deriving (Eq, Read, Show)++-- | O(/n/) Collect the /n/ simple powers of a sample.+--+-- Functions computed over a sample's simple powers require at least a+-- certain number (or /order/) of powers to be collected.+--+-- * To compute the /k/th 'centralMoment', at least /k/ simple powers+--   must be collected.+--+-- * For the 'variance', at least 2 simple powers are needed.+--+-- * For 'skewness', we need at least 3 simple powers.+--+-- * For 'kurtosis', at least 4 simple powers are required.+--+-- This function is subject to stream fusion.+powers :: Int                   -- ^ /n/, the number of powers, where /n/ >= 2.+       -> Sample+       -> Powers+powers k+    | k < 2     = error "Statistics.Sample.powers: too few powers"+    | otherwise = fini . foldlU go (unsafePerformIO . unsafeSTToIO $ create)+  where+    go ms x = inlinePerformIO . unsafeSTToIO $ loop 0 1+        where loop !i !xk | i == l = return ms+                          | otherwise = do+                readMU ms i >>= writeMU ms i . (+ xk)+                loop (i+1) (xk*x)+    fini = Powers . unsafePerformIO . unsafeSTToIO . unsafeFreezeAllMU+    create = newMU l >>= fill 0+        where fill !i ms | i == l    = return ms+                         | otherwise = writeMU ms i 0 >> fill (i+1) ms+    l = k + 1+{-# INLINE powers #-}++-- | The order (number) of simple powers collected from a 'Sample'.+order :: Powers -> Int+order (Powers pa) = lengthU pa - 1+{-# INLINE order #-}++-- | Compute the /k/th central moment of a 'Sample'.  The central+-- moment is also known as the moment about the mean.+centralMoment :: Int -> Powers -> Double+centralMoment k p@(Powers pa)+    | k < 0 || k > order p =+                  error ("Statistics.Sample.Powers.centralMoment: "+                         ++ "invalid argument")+    | k == 0    = 1+    | otherwise = (/n) . sumU . mapU go . indexedU . takeU (k+1) $ pa+  where+    go (i :*: e) = fromIntegral (k `choose` i) * ((-m) ^ (k-i)) * e+    n = indexU pa 0+    m = mean p+{-# INLINE centralMoment #-}++-- | Maximum likelihood estimate of a sample's variance.  Also known+-- as the population variance, where the denominator is /n/.  This is+-- the second central moment of the sample.+--+-- This is less numerically robust than the variance function in the+-- 'Statistics.Sample' module, but the number is essentially free to+-- compute if you have already collected a sample's simple powers.+--+-- Requires 'Powers' with 'order' at least 2.+variance :: Powers -> Double+variance = centralMoment 2+{-# INLINE variance #-}++-- | Standard deviation.  This is simply the square root of the+-- maximum likelihood estimate of the variance.+stdDev :: Powers -> Double+stdDev = sqrt . variance+{-# INLINE stdDev #-}++-- | Unbiased estimate of a sample's variance.  Also known as the+-- sample variance, where the denominator is /n/-1.+--+-- Requires 'Powers' with 'order' at least 2.+varianceUnbiased :: Powers -> Double+varianceUnbiased p@(Powers pa)+    | n > 1     = variance p * n / (n-1)+    | otherwise = 0+  where n = indexU pa 0+{-# INLINE varianceUnbiased #-}++-- | Compute the skewness of a sample. This is a measure of the+-- asymmetry of its distribution.+--+-- A sample with negative skew is said to be /left-skewed/.  Most of+-- its mass is on the right of the distribution, with the tail on the+-- left.+--+-- > skewness . powers 3 $ toU [1,100,101,102,103]+-- > ==> -1.497681449918257+--+-- A sample with positive skew is said to be /right-skewed/.+--+-- > skewness . powers 3 $ toU [1,2,3,4,100]+-- > ==> 1.4975367033335198+--+-- A sample's skewness is not defined if its 'variance' is zero.+--+-- Requires 'Powers' with 'order' at least 3.+skewness :: Powers -> Double+skewness p = centralMoment 3 p * variance p ** (-1.5)+{-# INLINE skewness #-}++-- | Compute the excess kurtosis of a sample.  This is a measure of+-- the \"peakedness\" of its distribution.  A high kurtosis indicates+-- that the sample's variance is due more to infrequent severe+-- deviations than to frequent modest deviations.+--+-- A sample's excess kurtosis is not defined if its 'variance' is+-- zero.+--+-- Requires 'Powers' with 'order' at least 4.+kurtosis :: Powers -> Double+kurtosis p = centralMoment 4 p / (v * v) - 3+    where v = variance p+{-# INLINE kurtosis #-}++-- | The number of elements in the original 'Sample'.  This is the+-- sample's zeroth simple power.+count :: Powers -> Int+count (Powers pa) = floor $ indexU pa 0+{-# INLINE count #-}++-- | The sum of elements in the original 'Sample'.  This is the+-- sample's first simple power.+sum :: Powers -> Double+sum (Powers pa) = indexU pa 1+{-# INLINE sum #-}++-- | The arithmetic mean of elements in the original 'Sample'.+--+-- This is less numerically robust than the mean function in the+-- 'Statistics.Sample' module, but the number is essentially free to+-- compute if you have already collected a sample's simple powers.+mean :: Powers -> Double+mean p@(Powers pa)+    | n == 0    = 0+    | otherwise = sum p / n+    where n     = indexU pa 0+{-# INLINE mean #-}++-- $references+--+-- * Besset, D.H. (2000) Elements of statistics.+--   /Object-oriented implementation of numerical methods/+--   pp. 311&#8211;331. <http://www.elsevier.com/wps/product/cws_home/677916>+--+-- * Anderson, G. (2009) Compute /k/th central moments in one+--   pass. /quantblog/. <http://quantblog.wordpress.com/2009/02/07/compute-kth-central-moments-in-one-pass/>

statistics.cabal view

@@ -1,5 +1,5 @@ name:           statistics-version:        0.2.1+version:        0.2.2 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful@@ -45,7 +45,10 @@     Statistics.Resampling     Statistics.Resampling.Bootstrap     Statistics.Sample+    Statistics.Sample.Powers     Statistics.Types+  other-modules:+    Statistics.Internal   build-depends:     base < 5,     erf,