statistics 0.5.1.2 → 0.6.0.0
raw patch · 9 files changed
+151/−127 lines, 9 files
Files
- Statistics/Autocorrelation.hs +13/−12
- Statistics/Function.hs +18/−17
- Statistics/KernelDensity.hs +23/−16
- Statistics/Math.hs +15/−11
- Statistics/Quantile.hs +23/−21
- Statistics/Resampling/Bootstrap.hs +1/−1
- Statistics/Sample.hs +46/−35
- Statistics/Sample/Powers.hs +11/−13
- statistics.cabal +1/−1
Statistics/Autocorrelation.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE FlexibleContexts #-} -- | -- Module : Statistics.Autocorrelation -- Copyright : (c) 2009 Bryan O'Sullivan@@ -16,31 +17,31 @@ , autocorrelation ) where -import Statistics.Sample (Sample, mean)-import qualified Data.Vector.Unboxed as U+import Statistics.Sample (mean)+import qualified Data.Vector.Generic as G -- | Compute the autocovariance of a sample, i.e. the covariance of -- the sample against a shifted version of itself.-autocovariance :: Sample -> U.Vector Double-autocovariance a = U.map f . U.enumFromTo 0 $ l-2+autocovariance :: (G.Vector v Double, G.Vector v Int) => v Double -> v Double+autocovariance a = G.map f . G.enumFromTo 0 $ l-2 where- f k = U.sum (U.zipWith (*) (U.take (l-k) c) (U.slice k (l-k) c))+ f k = G.sum (G.zipWith (*) (G.take (l-k) c) (G.slice k (l-k) c)) / fromIntegral l- c = U.map (subtract (mean a)) a- l = U.length a+ c = G.map (subtract (mean a)) a+ l = G.length 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 -> (U.Vector Double, U.Vector Double, U.Vector Double)+autocorrelation :: (G.Vector v Double, G.Vector v Int) => v Double -> (v Double, v Double, v Double) autocorrelation a = (r, ci (-), ci (+)) where- r = U.map (/ U.head c) c+ r = G.map (/ G.head c) c where c = autocovariance a- dllse = U.map f . U.scanl1 (+) . U.map square $ r+ dllse = G.map f . G.scanl1 (+) . G.map square $ r where f v = 1.96 * sqrt ((v * 2 + 1) / l)- l = fromIntegral (U.length a)- ci f = U.cons 1 . U.tail . U.map (f (-1/l)) $ dllse+ l = fromIntegral (G.length a)+ ci f = G.cons 1 . G.tail . G.map (f (-1/l)) $ dllse square x = x * x
Statistics/Function.hs view
@@ -1,4 +1,5 @@-{-# LANGUAGE Rank2Types, TypeOperators #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE FlexibleContexts #-} -- | -- Module : Statistics.Function -- Copyright : (c) 2009, 2010 Bryan O'Sullivan@@ -26,38 +27,38 @@ import Data.Vector.Algorithms.Combinators (apply) import Data.Vector.Generic (unsafeFreeze) import qualified Data.Vector.Algorithms.Intro as I-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Unboxed.Mutable as MU+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as M -- | Sort a vector.-sort :: (U.Unbox e, Ord e) => U.Vector e -> U.Vector e+sort :: (Ord e, G.Vector v e) => v e -> v e sort = apply I.sort {-# INLINE sort #-} -- | Partially sort a vector, such that the least /k/ elements will be -- at the front.-partialSort :: (U.Unbox e, Ord e) =>- Int -- ^ The number /k/ of least elements.- -> U.Vector e- -> U.Vector e+partialSort :: (G.Vector v e, Ord e) =>+ Int -- ^ The number /k/ of least elements.+ -> v e+ -> v e partialSort k = apply (\a -> I.partialSort a k) {-# INLINE partialSort #-} -- | Return the indices of a vector.-indices :: (U.Unbox a) => U.Vector a -> U.Vector Int-indices a = U.enumFromTo 0 (U.length a - 1)+indices :: (G.Vector v a, G.Vector v Int) => v a -> v Int+indices a = G.enumFromTo 0 (G.length a - 1) {-# INLINE indices #-} -- | Zip a vector with its indices.-indexed :: U.Unbox e => U.Vector e -> U.Vector (Int,e)-indexed a = U.zip (indices a) a+indexed :: (G.Vector v e, G.Vector v Int, G.Vector v (Int,e)) => v e -> v (Int,e)+indexed a = G.zip (indices a) a {-# INLINE indexed #-} data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double -- | Compute the minimum and maximum of a vector in one pass.-minMax :: U.Vector Double -> (Double , Double)-minMax = fini . U.foldl go (MM (1/0) (-1/0))+minMax :: (G.Vector v Double) => v Double -> (Double, Double)+minMax = fini . G.foldl' 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)@@ -65,12 +66,12 @@ -- | Create a vector, using the given action to populate each -- element.-create :: (PrimMonad m, U.Unbox e) => Int -> (Int -> m e) -> m (U.Vector e)+create :: (PrimMonad m, G.Vector v e) => Int -> (Int -> m e) -> m (v e) create size itemAt = assert (size >= 0) $- MU.new size >>= loop 0+ M.new size >>= loop 0 where loop k arr | k >= size = unsafeFreeze arr | otherwise = do r <- itemAt k- MU.write arr k r+ M.write arr k r loop (k+1) arr {-# INLINE create #-}
Statistics/KernelDensity.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE FlexibleContexts #-} -- | -- Module : Statistics.KernelDensity -- Copyright : (c) 2009 Bryan O'Sullivan@@ -40,8 +41,8 @@ import Statistics.Constants (m_1_sqrt_2, m_2_sqrt_pi) import Statistics.Function (minMax) import Statistics.Sample (stdDev)-import Statistics.Types (Sample) import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic as G -- | Points from the range of a 'Sample'. newtype Points = Points {@@ -61,10 +62,11 @@ -- | Compute the optimal bandwidth from the observed data for the given -- kernel.-bandwidth :: (Double -> Bandwidth)- -> Sample+bandwidth :: G.Vector v Double =>+ (Double -> Bandwidth)+ -> v Double -> Bandwidth-bandwidth kern values = stdDev values * kern (fromIntegral $ U.length values)+bandwidth kern values = stdDev values * kern (fromIntegral $ G.length values) -- | Choose a uniform range of points at which to estimate a sample's -- probability density function.@@ -74,9 +76,10 @@ -- -- If this function is passed an empty vector, it returns values of -- positive and negative infinity.-choosePoints :: Int -- ^ Number of points to select, /n/+choosePoints :: G.Vector v Double =>+ Int -- ^ Number of points to select, /n/ -> Double -- ^ Sample bandwidth, /h/- -> Sample -- ^ Input data+ -> v Double -- ^ Input data -> Points choosePoints n h sample = Points . U.map f $ U.enumFromTo 0 n' where lo = a - h@@ -116,27 +119,29 @@ -- | Kernel density estimator, providing a non-parametric way of -- estimating the PDF of a random variable.-estimatePDF :: Kernel -- ^ Kernel function+estimatePDF :: G.Vector v Double =>+ Kernel -- ^ Kernel function -> Bandwidth -- ^ Bandwidth, /h/- -> Sample -- ^ Sample data+ -> v Double -- ^ Sample data -> Points -- ^ Points at which to estimate -> U.Vector Double estimatePDF kernel h sample | n < 2 = errorShort "estimatePDF" | otherwise = U.map k . fromPoints where- k p = U.sum . U.map (kernel f h p) $ sample+ k p = G.sum . G.map (kernel f h p) $ sample f = 1 / (h * fromIntegral n)- n = U.length sample+ n = G.length 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+simplePDF :: G.Vector v Double =>+ (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+ -> v Double -- ^ sample data -> (Points, U.Vector Double) simplePDF fbw fpdf k numPoints sample = (points, estimatePDF fpdf bw sample points)@@ -147,16 +152,18 @@ -- | 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+epanechnikovPDF :: G.Vector v Double =>+ Int -- ^ Number of points at which to estimate+ -> v Double -- ^ Data sample -> (Points, U.Vector 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+gaussianPDF :: G.Vector v Double =>+ Int -- ^ Number of points at which to estimate+ -> v Double -- ^ Data sample -> (Points, U.Vector Double) gaussianPDF = simplePDF gaussianBW gaussianKernel 3
Statistics/Math.hs view
@@ -1,4 +1,5 @@-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-} -- | -- Module : Statistics.Math -- Copyright : (c) 2009 Bryan O'Sullivan@@ -26,25 +27,28 @@ -- $references ) where -import Data.Vector.Unboxed ((!))+import Data.Vector.Generic ((!)) import Data.Word (Word64) import Statistics.Constants (m_sqrt_2_pi) import Statistics.Distribution (cumulative) import Statistics.Distribution.Normal (standard) import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic as G data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double -- | Evaluate a series of Chebyshev polynomials. Uses Clenshaw's -- algorithm.-chebyshev :: Double -- ^ Parameter of each function.- -> U.Vector Double -- ^ Coefficients of each polynomial- -- term, in increasing order.+chebyshev :: (G.Vector v Double) =>+ Double -- ^ Parameter of each function.+ -> v Double -- ^ Coefficients of each polynomial term, in increasing order. -> Double-chebyshev x a = fini . U.foldl step (C 0 0) $ U.enumFromStepN (U.length a - 1) (-1) (U.length a - 1)+chebyshev x a = fini . U.foldl' step (C 0 0) $ U.enumFromStepN (len - 1) (-1) (len - 1) where step (C b1 b2) k = C ((a ! k) + x2 * b1 - b2) b1 fini (C b1 b2) = (a ! 0) + x * b1 - b2- x2 = x * 2+ x2 = x * 2+ len = G.length a+{-# INLINE chebyshev #-} -- | The binomial coefficient. --@@ -52,7 +56,7 @@ choose :: Int -> Int -> Double n `choose` k | k > n = 0- | k < 30 = U.foldl go 1 . U.enumFromTo 1 $ k'+ | k < 30 = U.foldl' go 1 . U.enumFromTo 1 $ k' | otherwise = exp $ lg (n+1) - lg (k+1) - lg (n-k+1) where go a i = a * (nk + j) / j where j = fromIntegral i :: Double@@ -71,8 +75,8 @@ factorial n | n < 0 = error "Statistics.Math.factorial: negative input" | n <= 1 = 0- | n <= 14 = fini . U.foldl goLong (F 1 1) $ ns- | otherwise = U.foldl goDouble 1 $ ns+ | n <= 14 = fini . U.foldl' goLong (F 1 1) $ ns+ | otherwise = U.foldl' goDouble 1 $ ns where goDouble t k = t * fromIntegral k goLong (F z x) _ = F (z * x') x' where x' = x + 1@@ -204,7 +208,7 @@ logGammaL :: Double -> Double logGammaL x | x <= 0 = 1/0- | otherwise = fini . U.foldl go (L 0 (x+7)) $ a+ | otherwise = fini . U.foldl' 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
Statistics/Quantile.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-} -- | -- Module : Statistics.Quantile -- Copyright : (c) 2009 Bryan O'Sullivan@@ -38,27 +38,27 @@ ) where import Control.Exception (assert)-import Data.Vector.Unboxed ((!))+import Data.Vector.Generic ((!)) import Statistics.Constants (m_epsilon) import Statistics.Function (partialSort)-import Statistics.Types (Sample)-import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic as G -- | 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.+weightedAvg :: G.Vector v Double => + Int -- ^ /k/, the desired quantile.+ -> Int -- ^ /q/, the number of quantiles.+ -> v Double -- ^ /x/, the sample data. -> Double weightedAvg k q x = assert (q >= 2) . assert (k >= 0) . assert (k < q) .- assert (U.all (not . isNaN) x) $+ assert (G.all (not . isNaN) x) $ xj + g * (xj1 - xj) where j = floor idx- idx = fromIntegral (U.length x - 1) * fromIntegral k / fromIntegral q+ idx = fromIntegral (G.length x - 1) * fromIntegral k / fromIntegral q g = idx - fromIntegral j xj = sx ! j xj1 = sx ! (j+1)@@ -72,16 +72,17 @@ -- 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.+continuousBy :: G.Vector v Double =>+ ContParam -- ^ Parameters /a/ and /b/.+ -> Int -- ^ /k/, the desired quantile.+ -> Int -- ^ /q/, the number of quantiles.+ -> v Double -- ^ /x/, the sample data. -> Double continuousBy (ContParam a b) k q x = assert (q >= 2) . assert (k >= 0) . assert (k <= q) .- assert (U.all (not . isNaN) x) $+ assert (G.all (not . isNaN) x) $ (1-h) * item (j-1) + h * item j where j = floor (t + eps)@@ -91,7 +92,7 @@ | otherwise = r where r = t - fromIntegral j eps = m_epsilon * 4- n = U.length x+ n = G.length x item = (sx !) . bracket sx = partialSort (bracket j + 1) x bracket m = min (max m 0) (n - 1)@@ -104,14 +105,15 @@ -- For instance, the interquartile range (IQR) can be estimated as -- follows: ----- > midspread medianUnbiased 4 (U.to [1,1,2,2,3])+-- > midspread medianUnbiased 4 (U.fromList [1,1,2,2,3]) -- > ==> 1.333333-midspread :: ContParam -- ^ Parameters /a/ and /b/.- -> Int -- ^ /q/, the number of quantiles.- -> Sample -- ^ /x/, the sample data.+midspread :: G.Vector v Double =>+ ContParam -- ^ Parameters /a/ and /b/.+ -> Int -- ^ /q/, the number of quantiles.+ -> v Double -- ^ /x/, the sample data. -> Double midspread (ContParam a b) k x =- assert (U.all (not . isNaN) x) .+ assert (G.all (not . isNaN) x) . assert (k > 0) $ quantile (1-frac) - quantile frac where@@ -122,7 +124,7 @@ | otherwise = r where r = t i - fromIntegral (j i) eps = m_epsilon * 4- n = U.length x+ n = G.length x item = (sx !) . bracket sx = partialSort (bracket (j (1-frac)) + 1) x bracket m = min (max m 0) (n - 1)
Statistics/Resampling/Bootstrap.hs view
@@ -83,7 +83,7 @@ ni = U.length resample n = fromIntegral ni accel = sumCubes / (6 * (sumSquares ** 1.5))- where (sumSquares :< sumCubes) = U.foldl f (0 :< 0) jack+ where (sumSquares :< sumCubes) = U.foldl' f (0 :< 0) jack f (s :< c) j = s + d2 :< c + d2 * d where d = jackMean - j d2 = d * d
Statistics/Sample.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-} -- | -- Module : Statistics.Sample -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan@@ -15,6 +15,7 @@ ( -- * Types Sample+ , WeightedSample -- * Descriptive functions , range @@ -52,18 +53,20 @@ import Statistics.Function (minMax) import Statistics.Types (Sample,WeightedSample)-import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic as G +-- Operator ^ will be overriden+import Prelude hiding ((^)) -range :: Sample -> Double+range :: (G.Vector v Double) => v Double -> Double range s = hi - lo where (lo , hi) = 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-mean = fini . U.foldl go (T 0 0)+mean :: (G.Vector v Double) => v Double -> Double+mean = fini . G.foldl' go (T 0 0) where fini (T a _) = a go (T m n) x = T m' n'@@ -73,8 +76,8 @@ -- | Arithmetic mean for weighted sample. It uses algorithm analogous -- to one in 'mean'-meanWeighted :: WeightedSample -> Double-meanWeighted = fini . U.foldl go (V 0 0)+meanWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double+meanWeighted = fini . G.foldl' go (V 0 0) where fini (V a _) = a go (V m w) (x,xw) = V m' w'@@ -85,16 +88,16 @@ -- | Harmonic mean. This algorithm performs a single pass over the -- sample.-harmonicMean :: Sample -> Double-harmonicMean = fini . U.foldl go (T 0 0)+harmonicMean :: (G.Vector v Double) => v Double -> Double+harmonicMean = fini . G.foldl' go (T 0 0) where 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 = fini . U.foldl go (T 1 0)+geometricMean :: (G.Vector v Double) => v Double -> Double+geometricMean = fini . G.foldl' go (T 1 0) where fini (T p n) = p ** (1 / fromIntegral n) go (T p n) a = T (p * a) (n + 1)@@ -108,12 +111,12 @@ -- -- For samples containing many values very close to the mean, this -- function is subject to inaccuracy due to catastrophic cancellation.-centralMoment :: Int -> Sample -> Double+centralMoment :: (G.Vector v Double) => Int -> v Double -> Double centralMoment a xs | a < 0 = error "Statistics.Sample.centralMoment: negative input" | a == 0 = 1 | a == 1 = 0- | otherwise = U.sum (U.map go xs) / fromIntegral (U.length xs)+ | otherwise = G.sum (G.map go xs) / fromIntegral (G.length xs) where go x = (x-m) ^ a m = mean xs@@ -126,15 +129,15 @@ -- -- 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 :: (G.Vector v Double) => Int -> Int -> v Double -> (Double, Double) centralMoments a b xs | a < 2 || b < 2 = (centralMoment a xs , centralMoment b xs)- | otherwise = fini . U.foldl go (V 0 0) $ xs+ | otherwise = fini . G.foldl' 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 (U.length xs)+ n = fromIntegral (G.length xs) {-# INLINE centralMoments #-} -- | Compute the skewness of a sample. This is a measure of the@@ -159,7 +162,7 @@ -- -- For samples containing many values very close to the mean, this -- function is subject to inaccuracy due to catastrophic cancellation.-skewness :: Sample -> Double+skewness :: (G.Vector v Double) => v Double -> Double skewness xs = c3 * c2 ** (-1.5) where (c3 , c2) = centralMoments 3 2 xs {-# INLINE skewness #-}@@ -177,7 +180,7 @@ -- -- For samples containing many values very close to the mean, this -- function is subject to inaccuracy due to catastrophic cancellation.-kurtosis :: Sample -> Double+kurtosis :: (G.Vector v Double) => v Double -> Double kurtosis xs = c4 / (c2 * c2) - 3 where (c4 , c2) = centralMoments 4 2 xs {-# INLINE kurtosis #-}@@ -201,48 +204,50 @@ sqr :: Double -> Double sqr x = x * x -robustSumVar :: Sample -> Double-robustSumVar samp = U.sum . U.map (sqr . subtract m) $ samp+robustSumVar :: (G.Vector v Double) => v Double -> Double+robustSumVar samp = G.sum . G.map (sqr . subtract m) $ samp where m = mean samp+{-# INLINE robustSumVar #-} -- | Maximum likelihood estimate of a sample's variance. Also known -- as the population variance, where the denominator is /n/.-variance :: Sample -> Double+variance :: (G.Vector v Double) => v Double -> Double variance samp | n > 1 = robustSumVar samp / fromIntegral n | otherwise = 0 where- n = U.length samp+ n = G.length samp {-# INLINE variance #-} -- | Unbiased estimate of a sample's variance. Also known as the -- sample variance, where the denominator is /n/-1.-varianceUnbiased :: Sample -> Double+varianceUnbiased :: (G.Vector v Double) => v Double -> Double varianceUnbiased samp | n > 1 = robustSumVar samp / fromIntegral (n-1) | otherwise = 0 where- n = U.length samp+ n = G.length samp {-# INLINE varianceUnbiased #-} -- | Standard deviation. This is simply the square root of the -- unbiased estimate of the variance.-stdDev :: Sample -> Double+stdDev :: (G.Vector v Double) => v Double -> Double stdDev = sqrt . varianceUnbiased-+{-# INLINE stdDev #-} -robustSumVarWeighted :: WeightedSample -> V-robustSumVarWeighted samp = U.foldl go (V 0 0) samp+robustSumVarWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> V+robustSumVarWeighted samp = G.foldl' go (V 0 0) samp where go (V s w) (x,xw) = V (s + xw*d*d) (w + xw) where d = x - m m = meanWeighted samp+{-# INLINE robustSumVarWeighted #-} -- | Weighted variance. This is biased estimation.-varianceWeighted :: WeightedSample -> Double+varianceWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double varianceWeighted samp- | U.length samp > 1 = fini $ robustSumVarWeighted samp+ | G.length samp > 1 = fini $ robustSumVarWeighted samp | otherwise = 0 where fini (V s w) = s / w@@ -259,8 +264,8 @@ -- mean, Knuth's algorithm gives inaccurate results due to -- catastrophic cancellation. -fastVar :: Sample -> T1-fastVar = U.foldl go (T1 0 0 0)+fastVar :: (G.Vector v Double) => v Double -> T1+fastVar = G.foldl' go (T1 0 0 0) where go (T1 n m s) x = T1 n' m' s' where n' = n + 1@@ -269,7 +274,7 @@ d = x - m -- | Maximum likelihood estimate of a sample's variance.-fastVariance :: Sample -> Double+fastVariance :: (G.Vector v Double) => v Double -> Double fastVariance = fini . fastVar where fini (T1 n _m s) | n > 1 = s / fromIntegral n@@ -277,7 +282,7 @@ {-# INLINE fastVariance #-} -- | Unbiased estimate of a sample's variance.-fastVarianceUnbiased :: Sample -> Double+fastVarianceUnbiased :: (G.Vector v Double) => v Double -> Double fastVarianceUnbiased = fini . fastVar where fini (T1 n _m s) | n > 1 = s / fromIntegral (n - 1)@@ -286,12 +291,18 @@ -- | Standard deviation. This is simply the square root of the -- maximum likelihood estimate of the variance.-fastStdDev :: Sample -> Double+fastStdDev :: (G.Vector v Double) => v Double -> Double fastStdDev = sqrt . fastVariance {-# INLINE fastStdDev #-} ------------------------------------------------------------------------ -- Helper code. Monomorphic unpacked accumulators.++-- (^) operator from Prelude is just slow.+(^) :: Double -> Int -> Double+x ^ 1 = x+x ^ n = x * (x ^ (n-1))+{-# INLINE (^) #-} -- don't support polymorphism, as we can't get unboxed returns if we use it. data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int
Statistics/Sample/Powers.hs view
@@ -1,4 +1,5 @@-{-# LANGUAGE BangPatterns, TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE BangPatterns #-} -- | -- Module : Statistics.Sample.Powers -- Copyright : (c) 2009, 2010 Bryan O'Sullivan@@ -19,8 +20,7 @@ module Statistics.Sample.Powers ( -- * Types- Sample- , Powers+ Powers -- * Constructor , powers@@ -53,9 +53,9 @@ import Statistics.Function (indexed) import Statistics.Internal (inlinePerformIO) import Statistics.Math (choose)-import Statistics.Types (Sample) import System.IO.Unsafe (unsafePerformIO) import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed.Mutable as MU newtype Powers = Powers (U.Vector Double)@@ -76,12 +76,13 @@ -- * 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 :: G.Vector v Double =>+ Int -- ^ /n/, the number of powers, where /n/ >= 2.+ -> v Double -> Powers powers k | k < 2 = error "Statistics.Sample.powers: too few powers"- | otherwise = fini . U.foldl go (unsafePerformIO $ create)+ | otherwise = fini . G.foldl' go (unsafePerformIO $ MU.newWith l 0) where go ms x = inlinePerformIO $ loop 0 1 where loop !i !xk | i == l = return ms@@ -89,18 +90,15 @@ MU.read ms i >>= MU.write ms i . (+ xk) loop (i+1) (xk*x) fini = Powers . unsafePerformIO . unsafeFreeze- create = MU.new l >>= fill 0- where fill !i ms | i == l = return ms- | otherwise = MU.write ms i 0 >> fill (i+1) ms- l = k + 1+ l = k + 1 {-# INLINE powers #-} --- | The order (number) of simple powers collected from a 'Sample'.+-- | The order (number) of simple powers collected from a 'sample'. order :: Powers -> Int order (Powers pa) = U.length pa - 1 {-# INLINE order #-} --- | Compute the /k/th central moment of a 'Sample'. The central+-- | 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)
statistics.cabal view
@@ -1,5 +1,5 @@ name: statistics-version: 0.5.1.2+version: 0.6.0.0 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful