diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs --- a/Statistics/Distribution/Exponential.hs +++ b/Statistics/Distribution/Exponential.hs @@ -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 -- ^ λ (scale) parameter. diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs --- a/Statistics/Distribution/Geometric.hs +++ b/Statistics/Distribution/Geometric.hs @@ -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 ( diff --git a/Statistics/Internal.hs b/Statistics/Internal.hs new file mode 100644 --- /dev/null +++ b/Statistics/Internal.hs @@ -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 diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs --- a/Statistics/Quantile.hs +++ b/Statistics/Quantile.hs @@ -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 diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs --- a/Statistics/Sample.hs +++ b/Statistics/Sample.hs @@ -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—and hence the standard deviation—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) diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs new file mode 100644 --- /dev/null +++ b/Statistics/Sample/Powers.hs @@ -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–331. +-- +-- * Anderson, G. (2009) Compute /k/th central moments in one +-- pass. /quantblog/. diff --git a/statistics.cabal b/statistics.cabal --- a/statistics.cabal +++ b/statistics.cabal @@ -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,