foldl-statistics (empty) → 0.1.0.0
raw patch · 6 files changed
+842/−0 lines, 6 filesdep +basedep +criteriondep +foldlsetup-changed
Dependencies added: base, criterion, foldl, foldl-statistics, math-functions, mwc-random, profunctors, quickcheck-instances, statistics, tasty, tasty-quickcheck, vector
Files
- LICENSE +112/−0
- Setup.hs +2/−0
- bench/Main.hs +135/−0
- foldl-statistics.cabal +63/−0
- src/Control/Foldl/Statistics.hs +398/−0
- test/Spec.hs +132/−0
+ LICENSE view
@@ -0,0 +1,112 @@+Copyright (c) 2009, 2010 Bryan O'Sullivan+Copyright (c) 2016, Commonwealth Scientific and Industrial Research Organisation+(CSIRO) ABN 41 687 119 230.++All rights reserved. CSIRO is willing to grant you a license to this+aemo-webservice on the following terms, except where otherwise indicated for+third party material.++Redistribution and use of this software 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.++* Neither the name of CSIRO nor the names of its contributors may be used to+ endorse or promote products derived from this software without specific prior+ written permission of CSIRO.++EXCEPT AS EXPRESSLY STATED IN THIS AGREEMENT AND TO THE FULL EXTENT PERMITTED BY+APPLICABLE LAW, THE SOFTWARE IS PROVIDED "AS-IS". CSIRO MAKES NO+REPRESENTATIONS, WARRANTIES OR CONDITIONS OF ANY KIND, EXPRESS OR IMPLIED,+INCLUDING BUT NOT LIMITED TO ANY REPRESENTATIONS, WARRANTIES OR CONDITIONS+REGARDING THE CONTENTS OR ACCURACY OF THE SOFTWARE, OR OF TITLE,+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NON-INFRINGEMENT, THE ABSENCE+OF LATENT OR OTHER DEFECTS, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT+DISCOVERABLE.++TO THE FULL EXTENT PERMITTED BY APPLICABLE LAW, IN NO EVENT SHALL CSIRO BE+LIABLE ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION, IN AN ACTION FOR+BREACH OF CONTRACT, NEGLIGENCE OR OTHERWISE) FOR ANY CLAIM, LOSS, DAMAGES OR+OTHER LIABILITY HOWSOEVER INCURRED. WITHOUT LIMITING THE SCOPE OF THE PREVIOUS+SENTENCE THE EXCLUSION OF LIABILITY SHALL INCLUDE: LOSS OF PRODUCTION OR+OPERATION TIME, LOSS, DAMAGE OR CORRUPTION OF DATA OR RECORDS; OR LOSS OF+ANTICIPATED SAVINGS, OPPORTUNITY, REVENUE, PROFIT OR GOODWILL, OR OTHER ECONOMIC+LOSS; OR ANY SPECIAL, INCIDENTAL, INDIRECT, CONSEQUENTIAL, PUNITIVE OR EXEMPLARY+DAMAGES, ARISING OUT OF OR IN CONNECTION WITH THIS AGREEMENT, ACCESS OF THE+SOFTWARE OR ANY OTHER DEALINGS WITH THE SOFTWARE, EVEN IF CSIRO HAS BEEN ADVISED+OF THE POSSIBILITY OF SUCH CLAIM, LOSS, DAMAGES OR OTHER LIABILITY.++APPLICABLE LEGISLATION SUCH AS THE AUSTRALIAN CONSUMER LAW MAY APPLY+REPRESENTATIONS, WARRANTIES, OR CONDITIONS, OR IMPOSES OBLIGATIONS OR LIABILITY+ON CSIRO THAT CANNOT BE EXCLUDED, RESTRICTED OR MODIFIED TO THE FULL EXTENT SET+OUT IN THE EXPRESS TERMS OF THIS CLAUSE ABOVE "CONSUMER GUARANTEES". TO THE+EXTENT THAT SUCH CONSUMER GUARANTEES CONTINUE TO APPLY, THEN TO THE FULL EXTENT+PERMITTED BY THE APPLICABLE LEGISLATION, THE LIABILITY OF CSIRO UNDER THE+RELEVANT CONSUMER GUARANTEE IS LIMITED (WHERE PERMITTED AT CSIRO'S OPTION) TO+ONE OF FOLLOWING REMEDIES OR SUBSTANTIALLY EQUIVALENT REMEDIES:++(a) THE REPLACEMENT OF THE SOFTWARE, THE SUPPLY OF EQUIVALENT+ SOFTWARE, OR SUPPLYING RELEVANT SERVICES AGAIN;+(b) THE REPAIR OF THE SOFTWARE;+(c) THE PAYMENT OF THE COST OF REPLACING THE+ SOFTWARE, OF ACQUIRING EQUIVALENT SOFTWARE, HAVING THE+ RELEVANT SERVICES SUPPLIED AGAIN, OR HAVING THE SOFTWARE+ REPAIRED.++IN THIS CLAUSE, CSIRO INCLUDES ANY THIRD PARTY AUTHOR OR OWNER OF ANY PART OF+THE SOFTWARE OR MATERIAL DISTRIBUTED WITH IT. CSIRO MAY ENFORCE ANY RIGHTS ON+BEHALF OF THE RELEVANT THIRD PARTY.++Third Party Components++The following third party components are distributed with the Software. You+agree to comply with the license terms for these components as part of accessing+the Software. Other third party software may also be identified in separate+files distributed with the Software.++___________________________________________________________________+___________________________________________________________________+The following Haskell library dependencies may be obtained from+https://hackage.haskell.org/packages/++StateVar 1.1.0.4+array 0.5.1.0+base 4.8.2.0+base-orphans 0.5.4+bifunctors 5.2+binary 0.7.5.0+bytestring 0.10.6.0+comonad 4.2.7.2+containers 0.5.6.2+contravariant 1.4+deepseq 1.4.1.1+distributive 0.5.0.2+erf 2.0.0.0+foldl 1.2.1+foldl-statistics 0.1.0.0+ghc-prim 0.4.0.0+hashable 1.2.4.0+integer-gmp 1.0.0.0+math-functions 0.1.7.0+mwc-random 0.13.4.0+primitive 0.6.1.0+profunctors 5.2+semigroups 0.18.1+stm 2.4.4.1+tagged 0.8.4+template-haskell 2.10.0.0+text 1.2.2.1+time 1.5.0.1+transformers 0.4.2.0+transformers-compat 0.4.0.4+unordered-containers 0.2.7.1+vector 0.11.0.0+vector-th-unbox 0.2.1.6+void 0.7.1+___________________________________________________________________
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bench/Main.hs view
@@ -0,0 +1,135 @@+module Main where++import Control.Monad.ST (runST)+import Criterion.Main+import qualified Statistics.Sample as S+import Statistics.Transform+import System.Random.MWC+import qualified Data.Vector.Unboxed as U+import Control.Foldl as F++import Control.Foldl.Statistics+++-- Test sample+sample :: U.Vector Double+sample = runST $ flip uniformVector 10000 =<< create++absSample = U.map abs sample++-- Weighted test sample+sampleW :: U.Vector (Double,Double)+sampleW = U.zip sample (U.reverse sample)++m = F.fold mean (U.toList sample)++mw = F.fold meanWeighted (U.toList sampleW)+++++main :: IO ()+main = defaultMain+ [ bgroup "Statistics of location"+ [ bgroup "mean"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold mean (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.mean sample+ ]+ , bgroup "meanWeighted"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold meanWeighted (U.toList vec)) sampleW+ , bench "Statistics.Sample" $ nf S.meanWeighted sampleW+ ]+ , bgroup "welfordMean"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold welfordMean (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.welfordMean sample+ ]+ , bgroup "harmonicMean"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold harmonicMean (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.harmonicMean sample+ ]+ , bgroup "geometricMean"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold geometricMean (U.toList vec)) absSample+ , bench "Statistics.Sample" $ nf S.geometricMean absSample+ ]+ ]+ , bgroup "Single-pass functions"+ [ bgroup "fastVariance"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold fastVariance (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.fastVariance sample+ ]+ , bgroup "fastVarianceUnbiased"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold fastVarianceUnbiased (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.fastVarianceUnbiased sample+ ]+ , bgroup "fastStdDev"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold fastStdDev (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.fastStdDev sample+ ]+ ]++ , bgroup "Functions requiring the mean"+ [ bgroup "variance"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (variance m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (variance (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.variance sample+ ]+ , bgroup "varianceUnbiased"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (varianceUnbiased m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (varianceUnbiased (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.varianceUnbiased sample+ ]+ , bgroup "stdDev"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (stdDev m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (stdDev (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.stdDev sample+ ]+ , bgroup "varianceWeighted"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (varianceWeighted m) (U.toList vec)) sampleW+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (varianceWeighted (F.fold meanWeighted (U.toList vec))) (U.toList vec)) sampleW+ , bench "Statistics.Sample" $ nf S.varianceWeighted sampleW+ ]+ ]++ , bgroup "Functions over central moments"+ [ bgroup "skewness"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (skewness m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (skewness (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.skewness sample+ ]+ , bgroup "kurtosis"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (kurtosis m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (kurtosis (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf S.kurtosis sample+ ]+ , bgroup "centralMoment 2"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (centralMoment 2 m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (centralMoment 2 (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf (S.centralMoment 2) sample+ ]+ , bgroup "centralMoment 3"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (centralMoment 3 m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (centralMoment 3 (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf (S.centralMoment 3) sample+ ]+ , bgroup "centralMoment 4"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (centralMoment 4 m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (centralMoment 4 (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf (S.centralMoment 4) sample+ ]+ , bgroup "centralMoment 7"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (centralMoment 7 m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (centralMoment 7 (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf (S.centralMoment 7) sample+ ]+ , bgroup "centralMoments 4 9"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (centralMoments 4 9 m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (centralMoments 4 9 (F.fold mean (U.toList vec))) (U.toList vec)) sample+ , bench "Statistics.Sample" $ nf (S.centralMoments 4 9) sample+ ]+ , bgroup "centralMoments' 4 9"+ [ bench "C.F.Statistics" $ nf (\vec -> F.fold (centralMoments' 4 9 m) (U.toList vec)) sample+ , bench "C.F.S(comp mean)" $ nf (\vec -> F.fold (centralMoments' 4 9 (F.fold mean (U.toList vec))) (U.toList vec)) sample+ ]+ ]+ ]+
+ foldl-statistics.cabal view
@@ -0,0 +1,63 @@+name: foldl-statistics+version: 0.1.0.0+synopsis: Statistical functions from the statistics package implemented as+ Folds.+description: The use of this package allows statistics to be computed using at most two+ passes over the input data, one to compute a mean and one to compute a further+ statistic such as variance and /n/th central moments. All algorithms are the+ obvious implementation of Bryan O\'Sullivan\'s+ <https://hackage.haskell.org/package/statistics statistics> package imeplemented+ as `Fold's from the+ <https://hackage.haskell.org/package/foldl foldl> package.+homepage: http://github.com/Data61/foldl-statistics#readme+license: BSD3+license-file: LICENSE+author: Alex Mason+maintainer: Alex.Mason@data61.csiro.au+copyright: 2016 Data61 (CSIRO)+category: Math, Statistics+build-type: Simple+-- extra-source-files:+cabal-version: >=1.10++library+ hs-source-dirs: src+ exposed-modules: Control.Foldl.Statistics+ default-language: Haskell2010+ build-depends: base >= 4.7 && < 5+ , foldl >= 1.1 && < 1.3+ , math-functions >= 0.1 && < 0.2+ , profunctors >= 5.2 && < 5.3++test-suite foldl-statistics-test+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Spec.hs+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010+ build-depends: base >= 4.7 && < 5.0+ , foldl-statistics+ , foldl >= 1.1 && < 1.3+ , statistics >= 0.13 && < 0.14+ , tasty >= 0.11 && < 0.12+ , tasty-quickcheck >= 0.8 && < 0.9+ , vector >= 0.11 && < 0.12+ , quickcheck-instances >= 0.3 && < 0.4+ , profunctors >= 5.2 && < 5.3++Benchmark bench-folds+ type: exitcode-stdio-1.0+ hs-source-dirs: bench+ main-is: Main.hs+ default-language: Haskell2010+ build-depends: base+ , foldl-statistics+ , criterion >= 1.1 && < 1.2+ , vector+ , statistics+ , mwc-random >= 0.13 && < 0.14+ , foldl++source-repository head+ type: git+ location: https://github.com/Data61/foldl-statistics
+ src/Control/Foldl/Statistics.hs view
@@ -0,0 +1,398 @@+-- |+-- Module : Control.Foldl.Statistics+-- Copyright : (c) 2011 Bryan O'Sullivan, 2016 National ICT Australia+-- License : BSD3+--+-- Maintainer : alex.mason@nicta.com.au+-- Stability : experimental+-- Portability : portable+--++module Control.Foldl.Statistics (+ -- * Introduction+ -- $intro+ -- * Descriptive functions+ range+ , sum'++ -- * Statistics of location+ , mean+ , welfordMean+ , meanWeighted+ , harmonicMean+ , geometricMean++ -- * Statistics of dispersion+ -- $variance++ -- ** Functions over central moments+ , centralMoment+ , centralMoments+ , centralMoments'+ , skewness+ , kurtosis++ -- ** Functions requiring the mean to be known (numerically robust)+ -- $robust+ , variance+ , varianceUnbiased+ , stdDev+ , varianceWeighted++ -- ** Single-pass functions (faster, less safe)+ -- $cancellation+ , fastVariance+ , fastVarianceUnbiased+ , fastStdDev++ -- $correlation+ , correlation++ -- * References+ -- $references+ , module Control.Foldl++ ) where++import Control.Foldl as F+import qualified Control.Foldl+import Data.Profunctor++import Numeric.Sum (KBNSum, kbn, add, zero)++data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int+data TS = TS {-# UNPACK #-}!KBNSum {-# UNPACK #-}!Int+data T1 = T1 {-# UNPACK #-}!Int {-# UNPACK #-}!Double {-# UNPACK #-}!Double+data V = V {-# UNPACK #-}!Double {-# UNPACK #-}!Double+data V1 = V1 {-# UNPACK #-}!Double {-# UNPACK #-}!Double {-# UNPACK #-}!Int+data V1S = V1S {-# UNPACK #-}!KBNSum {-# UNPACK #-}!KBNSum {-# UNPACK #-}!Int+++-- $intro+-- Statistical functions from the+-- <https://hackage.haskell.org/package/statistics/docs/Statistics-Sample.html Statistics.Sample>+-- module of the+-- <https://hackage.haskell.org/package/statistics statistics> package by+-- Bryan O'Sullivan, implemented as `Control.Foldl.Fold's from the+-- <https://hackage.haskell.org/package/foldl foldl> package.+--+-- This allows many statistics to be computed concurrently with at most+-- two passes over the data, usually by computing the `mean' first, and+-- passing it to further `Fold's.++-- | A numerically stable sum using Kahan-Babuška-Neumaier+-- summation from "Numeric.Sum"+{-# INLINE sum' #-}+sum' :: Fold Double Double+sum' = Fold (add :: KBNSum -> Double -> KBNSum)+ (zero :: KBNSum)+ kbn+++-- | The difference between the largest and smallest+-- elements of a sample.+{-# INLINE range #-}+range :: Fold Double Double+range = (\(Just lo) (Just hi) -> hi - lo)+ <$> F.minimum+ <*> F.maximum++-- | Arithmetic mean. This uses Kahan-Babuška-Neumaier+-- summation, so is more accurate than 'welfordMean' unless the input+-- values are very large.+{-# INLINE mean #-}+mean :: Fold Double Double+mean = Fold step (TS zero 0) final where+ step (TS s n) x = TS (add s x) (n+1)+ final (TS s n) = kbn s / fromIntegral n+++-- | Arithmetic mean. This uses Welford's algorithm to provide+-- numerical stability, using a single pass over the sample data.+--+-- Compared to 'mean', this loses a surprising amount of precision+-- unless the inputs are very large.+{-# INLINE welfordMean #-}+welfordMean :: Fold Double Double+welfordMean = Fold step (T 0 0) final where+ final (T m _) = m+ step (T m n) x = T m' n' where+ m' = m + (x - m) / fromIntegral n'+ n' = n + 1+++-- | Arithmetic mean for weighted sample. It uses a single-pass+-- algorithm analogous to the one used by 'welfordMean'.+{-# INLINE meanWeighted #-}+meanWeighted :: Fold (Double,Double) Double+meanWeighted = Fold step (V 0 0) final+ where+ final (V a _) = a+ step (V m w) (x,xw) = V m' w'+ where m' | w' == 0 = 0+ | otherwise = m + xw * (x - m) / w'+ w' = w + xw++-- | Harmonic mean.+{-# INLINE harmonicMean #-}+harmonicMean :: Fold Double Double+harmonicMean = Fold step (T 0 0) final+ where+ final (T b a) = fromIntegral a / b+ step (T x y) n = T (x + (1/n)) (y+1)++-- | Geometric mean of a sample containing no negative values.+{-# INLINE geometricMean #-}+geometricMean :: Fold Double Double+geometricMean = dimap log exp mean++-- | Compute the /k/th central moment of a sample. The central moment+-- is also known as the moment about the mean.+--+-- This function requires the mean of the data to compute the central moment.+--+-- For samples containing many values very close to the mean, this+-- function is subject to inaccuracy due to catastrophic cancellation.+{-# INLINE centralMoment #-}+centralMoment :: Int -> Double -> Fold Double Double+centralMoment a m+ | a < 0 = error "Statistics.Sample.centralMoment: negative input"+ | a == 0 = 1+ | a == 1 = 0+ | otherwise = Fold step (TS zero 0) final where+ step (TS s n) x = TS (add s $ go x) (n+1)+ final (TS s n) = kbn s / fromIntegral n+ go x = (x-m) ^^^ a++-- | Compute the /k/th and /j/th central moments of a sample.+--+-- This fold requires the mean of the data to be known.+--+-- For samples containing many values very close to the mean, this+-- function is subject to inaccuracy due to catastrophic cancellation.+{-# INLINE centralMoments #-}+centralMoments :: Int -> Int -> Double -> Fold Double (Double, Double)+centralMoments a b m+ | a < 2 || b < 2 = (,) <$> centralMoment a m <*> centralMoment b m+ | otherwise = Fold step (V1 0 0 0) final+ where final (V1 i j n) = (i / fromIntegral n , j / fromIntegral n)+ step (V1 i j n) x = V1 (i + d^^^a) (j + d^^^b) (n+1)+ where d = x - m+++-- | Compute the /k/th and /j/th central moments of a sample.+--+-- This fold requires the mean of the data to be known.+--+-- This variation of `centralMoments' uses Kahan-Babuška-Neumaier+-- summation to attempt to improve the accuracy of results, which may+-- make computation slower.+{-# INLINE centralMoments' #-}+centralMoments' :: Int -> Int -> Double -> Fold Double (Double, Double)+centralMoments' a b m+ | a < 2 || b < 2 = (,) <$> centralMoment a m <*> centralMoment b m+ | otherwise = Fold step (V1S zero zero 0) final+ where final (V1S i j n) = (kbn i / fromIntegral n , kbn j / fromIntegral n)+ step (V1S i j n) x = V1S (add i $ d^^^a) (add j $ d^^^b) (n+1)+ where d = x - m++-- | 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 $ U.to [1,100,101,102,103]+-- > ==> -1.497681449918257+--+-- A sample with positive skew is said to be /right-skewed/.+--+-- > skewness $ U.to [1,2,3,4,100]+-- > ==> 1.4975367033335198+--+-- A sample's skewness is not defined if its 'variance' is zero.+--+-- This fold requires the mean of the data to be known.+--+-- For samples containing many values very close to the mean, this+-- function is subject to inaccuracy due to catastrophic cancellation.+{-# INLINE skewness #-}+skewness :: Double -> Fold Double Double+skewness m = (\(c3, c2) -> c3 * c2 ** (-1.5)) <$> centralMoments 3 2 m+++-- | 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 fold requires the mean of the data to be known.+--+-- For samples containing many values very close to the mean, this+-- function is subject to inaccuracy due to catastrophic cancellation.+{-# INLINE kurtosis #-}+kurtosis :: Double -> Fold Double Double+kurtosis m = (\(c4,c2) -> c4 / (c2 * c2) - 3) <$> centralMoments 4 2 m+++-- $variance+--+-- The variance—and hence the standard deviation—of a+-- sample of fewer than two elements are both defined to be zero.+--+-- Many of these Folds take the mean as an argument for constructing+-- the variance, and as such require two passes over the data.++-- $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.+++-- Multiply a number by itself.+{-# INLINE square #-}+square :: Double -> Double+square x = x * x++{-# INLINE robustSumVar #-}+robustSumVar :: Double -> Fold Double TS+robustSumVar m = Fold step (TS zero 0) id where+ step (TS s n) x = TS (add s . square . subtract m $ x) (n+1)++-- | Maximum likelihood estimate of a sample's variance. Also known+-- as the population variance, where the denominator is /n/.+{-# INLINE variance #-}+variance :: Double -> Fold Double Double+variance m =+ (\(TS sv n) -> if n > 1 then kbn sv / fromIntegral n else 0)+ <$> robustSumVar m++-- | Unbiased estimate of a sample's variance. Also known as the+-- sample variance, where the denominator is /n/-1.+{-# INLINE varianceUnbiased #-}+varianceUnbiased :: Double -> Fold Double Double+varianceUnbiased m =+ (\(TS sv n) -> if n > 1 then kbn sv / fromIntegral (n-1) else 0)+ <$> robustSumVar m+++-- | Standard deviation. This is simply the square root of the+-- unbiased estimate of the variance.+{-# INLINE stdDev #-}+stdDev :: Double -> Fold Double Double+stdDev m = sqrt (varianceUnbiased m)+++{-# INLINE robustSumVarWeighted #-}+robustSumVarWeighted :: Double -> Fold (Double,Double) V1+robustSumVarWeighted m = Fold step (V1 0 0 0) id+ where+ step (V1 s w n) (x,xw) = V1 (s + xw*d*d) (w + xw) (n+1)+ where d = x - m++-- | Weighted variance. This is biased estimation. Requires the+-- weighted mean of the input data.+{-# INLINE varianceWeighted #-}+varianceWeighted :: Double -> Fold (Double,Double) Double+varianceWeighted m =+ (\(V1 s w n) -> if n > 1 then s / w else 0)+ <$> robustSumVarWeighted m++-- $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 fusion and do not require the mean to be passed.+--+-- /Note/: in cases where most sample data is close to the sample's+-- mean, Knuth's algorithm gives inaccurate results due to+-- catastrophic cancellation.++{-# INLINE fastVar #-}+fastVar :: Fold Double T1+fastVar = Fold step (T1 0 0 0) id+ where+ step (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.+{-# INLINE fastVariance #-}+fastVariance :: Fold Double Double+fastVariance = final <$> fastVar+ where final (T1 n _m s)+ | n > 1 = s / fromIntegral n+ | otherwise = 0+++-- | Maximum likelihood estimate of a sample's variance.+{-# INLINE fastVarianceUnbiased #-}+fastVarianceUnbiased :: Fold Double Double+fastVarianceUnbiased = final <$> fastVar+ where final (T1 n _m s)+ | n > 1 = s / fromIntegral (n-1)+ | otherwise = 0+++-- | Standard deviation. This is simply the square root of the+-- maximum likelihood estimate of the variance.+{-# INLINE fastStdDev #-}+fastStdDev :: Fold Double Double+fastStdDev = sqrt fastVariance+++-- $correlation+--+--+correlation :: (Double, Double) -> (Double, Double) -> Fold (Double,Double) Double+correlation (m1,m2) (s1,s2) = Fold step (TS zero 0) final where+ step (TS s n) (x1,x2) = TS (add s $ ((x1-m1)/s1) * ((x2-m2)/s2)) (n+1)+ final (TS s n) = kbn s / fromIntegral (n-1)+++-- $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–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–535. <http://doi.acm.org/10.1145/359146.359153>++++-- (^) operator from Prelude is just slow.+(^^^) :: Double -> Int -> Double+x ^^^ 1 = x+x ^^^ n = x * (x ^^^ (n-1))+{-# INLINE[2] (^^^) #-}+{-# RULES+"pow 2" forall x. x ^^^ 2 = x * x+"pow 3" forall x. x ^^^ 3 = x * x * x+"pow 4" forall x. x ^^^ 4 = x * x * x * x+"pow 5" forall x. x ^^^ 5 = x * x * x * x * x+"pow 6" forall x. x ^^^ 6 = x * x * x * x * x * x+"pow 7" forall x. x ^^^ 7 = x * x * x * x * x * x * x+"pow 8" forall x. x ^^^ 8 = x * x * x * x * x * x * x * x+"pow 9" forall x. x ^^^ 9 = x * x * x * x * x * x * x * x * x+"pow 10" forall x. x ^^^ 10 = x * x * x * x * x * x * x * x * x * x++ #-}
+ test/Spec.hs view
@@ -0,0 +1,132 @@++import Test.Tasty+-- import Test.Tasty.SmallCheck as SC+import qualified Test.Tasty.QuickCheck as QC+import Test.Tasty.QuickCheck ((==>))++import qualified Control.Foldl as F+import Control.Foldl.Statistics hiding (length)++import qualified Data.Vector.Unboxed as U+import Test.QuickCheck.Instances++import qualified Statistics.Sample as S+import Statistics.Function (within)++import Data.Profunctor+++toV :: [Double] -> U.Vector Double+toV = U.fromList+++onVec :: String -> (U.Vector Double -> QC.Property) -> TestTree+onVec str f = QC.testProperty str (f . toV)++onVec2 :: String -> (U.Vector (Double,Double) -> QC.Property) -> TestTree+onVec2 str f = QC.testProperty str (f . U.fromList)++main :: IO ()+main = defaultMain $+ testGroup "Results match Statistics.Sample"+ [ testGroup "Without pre-computed mean"+ [ testGroup "Statistics of location"+ [ onVec "mean" $ \vec ->+ not (U.null vec) ==> F.fold mean (U.toList vec) == S.mean vec+ , onVec2 "meanWeighted" $ \vec ->+ not (U.null vec) ==> F.fold meanWeighted (U.toList vec) == S.meanWeighted vec+ , onVec "welfordMean" $ \vec ->+ not (U.null vec) ==> F.fold welfordMean (U.toList vec) == S.welfordMean vec+ , onVec "harmonicMean" $ \vec ->+ not (U.null vec) ==> F.fold harmonicMean (U.toList vec) == S.harmonicMean vec+ , onVec "geometricMean" $ \vec ->+ not (U.null vec) ==> let vec' = U.map abs vec+ in F.fold geometricMean (U.toList vec') == S.geometricMean vec'+ ]++ , testGroup "Single-pass functions"+ [ onVec "fastVariance" $ \vec ->+ not (U.null vec) ==> F.fold fastVariance (U.toList vec) == S.fastVariance vec+ , onVec "fastVarianceUnbiased" $ \vec ->+ not (U.null vec) ==> F.fold fastVarianceUnbiased (U.toList vec) == S.fastVarianceUnbiased vec+ , onVec "fastStdDev" $ \vec ->+ not (U.null vec) ==> F.fold fastStdDev (U.toList vec) == S.fastStdDev vec+ ]+ ]++ , testGroup "With pre-computed mean"+ [ testGroup "Functions requiring the mean to be known"+ [ onVec "variance" $ \vec ->+ not (U.null vec) ==> let m = F.fold mean (U.toList vec)+ in F.fold (variance m) (U.toList vec) == S.variance vec+ , onVec "varianceUnbiased" $ \vec ->+ not (U.null vec) ==> let m = F.fold mean (U.toList vec)+ in F.fold (varianceUnbiased m) (U.toList vec) == S.varianceUnbiased vec+ , onVec "stdDev" $ \vec ->+ not (U.null vec) ==> let m = F.fold mean (U.toList vec)+ in F.fold (stdDev m) (U.toList vec) == S.stdDev vec+ , onVec2 "varianceWeighted" $ \vec ->+ not (U.null vec) ==> let m = F.fold meanWeighted (U.toList vec)+ in F.fold (varianceWeighted m) (U.toList vec) == S.varianceWeighted vec+ ]++ , testGroup "Functions over central moments"+ [ onVec "skewness" $ \vec ->+ U.length vec > 3 ==> let m = F.fold mean (U.toList vec)+ in F.fold (skewness m) (U.toList vec) == S.skewness vec+ , onVec "kurtosis" $ \vec ->+ U.length vec > 4 ==> let m = F.fold mean (U.toList vec)+ in F.fold (kurtosis m) (U.toList vec) == S.kurtosis vec+ , onVec "centralMoment 2" $ \vec ->+ U.length vec > 2 ==> let m = F.fold mean (U.toList vec)+ in F.fold (centralMoment 2 m) (U.toList vec) == S.centralMoment 2 vec+ , onVec "centralMoment 3" $ \vec ->+ U.length vec > 3 ==> let m = F.fold mean (U.toList vec)+ in F.fold (centralMoment 3 m) (U.toList vec) == S.centralMoment 3 vec+ , onVec "centralMoment 4" $ \vec ->+ U.length vec > 4 ==> let m = F.fold mean (U.toList vec)+ in F.fold (centralMoment 4 m) (U.toList vec) == S.centralMoment 4 vec+ , onVec "centralMoment 7" $ \vec ->+ U.length vec > 7 ==> let m = F.fold mean (U.toList vec)+ in F.fold (centralMoment 7 m) (U.toList vec) == S.centralMoment 7 vec+ , onVec "centralMoments 4 9" $ \vec ->+ U.length vec > 7 ==> let m = F.fold mean (U.toList vec)+ in F.fold (centralMoments 4 9 m) (U.toList vec) == S.centralMoments 4 9 vec+ -- Cannot test this because we do not have an equivalent implementation+ -- from the statistics package.+ -- , onVec "centralMoments' 4 9" $ \vec -> length lst > 7 ==>+ -- let m = F.fold mean lst+ -- (f1,f2) = (F.fold (centralMoments' 4 9 m) lst)+ -- (s1,s2) = (S.centralMoments 4 9 vec)+ -- in within 3 f1 s1 && within 3 f2 s2+ ]+ , testGroup "Correlation"+ [ onVec2 "correlation between [-1,1]" $ \vec ->+ U.length vec > 2 ==>+ let m1 = F.fold mean (U.toList $ U.map fst vec)+ m2 = F.fold mean (U.toList $ U.map snd vec)+ s1 = F.fold (stdDev m1) (U.toList $ U.map fst vec)+ s2 = F.fold (stdDev m2) (U.toList $ U.map snd vec)+ in between (-1,1) $+ F.fold (correlation (m1,m2) (s1,s2)) (U.toList vec)+ , onVec2 "correlation between [-1,1] fastStdDev" $ \vec ->++ let (m1,m2) = F.fold ((,)+ <$> lmap fst mean+ <*> lmap snd mean)+ (U.toList vec)+ (s1,s2) = F.fold ((,)+ <$> lmap fst (stdDev m1)+ <*> lmap snd (stdDev m2))+ (U.toList vec)+ corr = F.fold (correlation (m1,m2) (s1,s2)) (U.toList vec)+ in U.length vec > 2 && s2 /= 0.0 && s2 /= 0.0 ==>+ QC.counterexample ("Correlation: " ++ show corr ++ " Stats: " ++ show (m1,m2,s1,s2)) $+ between (-1,1) corr || isNaN corr++ ]+ ]+ ]++between :: (Double,Double) -> Double -> Bool+between (lo,hi) = \x -> lo <= x && x <= hi