diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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 +___________________________________________________________________ diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/bench/Main.hs b/bench/Main.hs new file mode 100644 --- /dev/null +++ b/bench/Main.hs @@ -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 + ] + ] + ] + diff --git a/foldl-statistics.cabal b/foldl-statistics.cabal new file mode 100644 --- /dev/null +++ b/foldl-statistics.cabal @@ -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 + package imeplemented + as `Fold's from the + 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 diff --git a/src/Control/Foldl/Statistics.hs b/src/Control/Foldl/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/Control/Foldl/Statistics.hs @@ -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 +-- +-- module of the +-- package by +-- Bryan O'Sullivan, implemented as `Control.Foldl.Fold's from the +-- 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. +-- +-- * 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. +-- +-- * West, D.H.D. (1979) Updating mean and variance estimates: an +-- improved method. /Communications of the ACM/ +-- 22(9):532–535. + + + +-- (^) 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 + + #-} diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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