diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,7 @@ # 0.1.4.0 +- foldl >= 1.2.2 exports `mean` and ` variance`, so hide them. + +# 0.1.4.0 - Added monoidal interface to linear regression # 0.1.3.0 diff --git a/bench/Main.hs b/bench/Main.hs --- a/bench/Main.hs +++ b/bench/Main.hs @@ -1,14 +1,20 @@ +{-# LANGUAGE CPP #-} + 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.Monad.ST (runST) +import Criterion.Main +import qualified Data.Vector.Unboxed as U +import qualified Statistics.Sample as S +-- import Statistics.Transform +import System.Random.MWC +#if MIN_VERSION_foldl(1,2,2) +import Control.Foldl as F hiding (mean, variance) +#else +import Control.Foldl as F +#endif -import Control.Foldl.Statistics +import Control.Foldl.Statistics -- Test sample @@ -21,6 +27,7 @@ sample2 = runST $ flip uniformVector 10000 =<< create {-# NOINLINE absSample #-} +absSample :: U.Vector Double absSample = U.map abs sample -- Weighted test sample @@ -28,6 +35,7 @@ sampleW :: U.Vector (Double,Double) sampleW = U.zip sample (U.reverse sample) +m, mw :: Double m = F.fold mean (U.toList sample) mw = F.fold meanWeighted (U.toList sampleW) diff --git a/foldl-statistics.cabal b/foldl-statistics.cabal --- a/foldl-statistics.cabal +++ b/foldl-statistics.cabal @@ -1,64 +1,72 @@ -name: foldl-statistics -version: 0.1.4.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: CHANGELOG.md, README.md -cabal-version: >=1.10 +name: foldl-statistics +version: 0.1.4.1 +cabal-version: >=1.10 +build-type: Simple +license: BSD3 +license-file: LICENSE +copyright: 2016 Data61 (CSIRO) +maintainer: Alex.Mason@data61.csiro.au +homepage: http://github.com/Data61/foldl-statistics#readme +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. +category: Math, Statistics +author: Alex Mason +extra-source-files: + CHANGELOG.md + README.md +source-repository head + type: git + location: https://github.com/Data61/foldl-statistics + 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.3 - , profunctors >= 5.2 && < 5.3 - , semigroups + exposed-modules: + Control.Foldl.Statistics + build-depends: + base >=4.7 && <5, + foldl >=1.1 && <1.3, + math-functions >=0.1 && <0.3, + profunctors ==5.2.*, + semigroups >=0.18.2 && <1.0 + default-language: Haskell2010 + hs-source-dirs: src 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 + type: exitcode-stdio-1.0 + main-is: Spec.hs + build-depends: + base >=4.7 && <5.0, + foldl-statistics >=0.1.4.1 && <0.2, + foldl >=1.2.1 && <1.3, + statistics ==0.13.*, + tasty ==0.11.*, + tasty-quickcheck ==0.8.*, + vector ==0.11.*, + quickcheck-instances ==0.3.*, + profunctors ==5.2.*, + semigroups >=0.18.2 && <0.19 + default-language: Haskell2010 + hs-source-dirs: test + ghc-options: -threaded -rtsopts -with-rtsopts=-N -source-repository head - type: git - location: https://github.com/Data61/foldl-statistics +benchmark bench-folds + type: exitcode-stdio-1.0 + main-is: Main.hs + build-depends: + base >=4.9.0.0 && <4.10, + foldl-statistics >=0.1.4.1 && <0.2, + criterion ==1.1.*, + vector >=0.10 && <1.0, + statistics >=0.13.3.0 && <0.14, + mwc-random ==0.13.*, + foldl >=1.2.1 && <1.3 + default-language: Haskell2010 + hs-source-dirs: bench diff --git a/src/Control/Foldl/Statistics.hs b/src/Control/Foldl/Statistics.hs --- a/src/Control/Foldl/Statistics.hs +++ b/src/Control/Foldl/Statistics.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} -- | -- Module : Control.Foldl.Statistics -- Copyright : (c) 2011 Bryan O'Sullivan, 2016 National ICT Australia @@ -58,6 +59,7 @@ , getLinRegResult , LinRegResult(..) , LinRegState + , lrrCount , correlation -- * References @@ -66,13 +68,22 @@ ) where -import Control.Foldl as F +#if MIN_VERSION_foldl(1,2,2) +import Control.Foldl as F hiding (mean, variance) +#else +import Control.Foldl as F +#endif + import qualified Control.Foldl -import Data.Profunctor -import Data.Semigroup +import Data.Profunctor +import Data.Semigroup -import Numeric.Sum (KBNSum, kbn, add, zero) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +#endif +import Numeric.Sum (KBNSum, add, kbn, zero) + data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int data TS = TS {-# UNPACK #-}!KBNSum {-# UNPACK #-}!Int data T1 = T1 {-# UNPACK #-}!Int {-# UNPACK #-}!Double {-# UNPACK #-}!Double @@ -113,6 +124,9 @@ -- | Arithmetic mean. This uses Kahan-Babuška-Neumaier -- summation, so is more accurate than 'welfordMean' unless the input -- values are very large. +-- +-- Since foldl-1.2.2, 'Control.Foldl` exports a `mean` function, so you will +-- have to hide one. {-# INLINE mean #-} mean :: Fold Double Double mean = Fold step (TS zero 0) final where @@ -175,7 +189,7 @@ | 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 + go x = (x - m) ^^^ a -- | Compute the /k/th and /j/th central moments of a sample. -- @@ -393,8 +407,8 @@ n = an+bn n2 = n*n nd = fi n - and = fi an - bnd = fi bn + nda = fi an + ndb = fi bn -- delta = b.M1 - a.M1; delta = bm1 - am1 -- delta2 = delta*delta; @@ -404,20 +418,20 @@ -- delta4 = delta2*delta2; delta4 = delta2*delta2 -- combined.M1 = (a.n*a.M1 + b.n*b.M1) / combined.n; - m1 = (and*am1 + bnd*bm1 ) / nd + m1 = (nda*am1 + ndb*bm1 ) / nd -- combined.M2 = a.M2 + b.M2 + delta2*a.n*b.n / combined.n; - m2 = am2 + bm2 + delta2*and*bnd / nd + m2 = am2 + bm2 + delta2*nda*ndb / nd -- combined.M3 = a.M3 + b.M3 + delta3*a.n*b.n* (a.n - b.n)/(combined.n*combined.n); - m3 = am3 + bm3 + delta3*and*bnd* fi( an - bn )/ fi n2 + m3 = am3 + bm3 + delta3*nda*ndb* fi( an - bn )/ fi n2 -- combined.M3 += 3.0*delta * (a.n*b.M2 - b.n*a.M2) / combined.n; - + 3.0*delta * (and*bm2 - bnd*am2 ) / nd + + 3.0*delta * (nda*bm2 - ndb*am2 ) / nd -- -- combined.M4 = a.M4 + b.M4 + delta4*a.n*b.n * (a.n*a.n - a.n*b.n + b.n*b.n) /(combined.n*combined.n*combined.n); - m4 = am4 + bm4 + delta4*and*bnd *fi(an*an - an*bn + bn*bn ) / fi (n*n*n) + m4 = am4 + bm4 + delta4*nda*ndb *fi(an*an - an*bn + bn*bn ) / fi (n*n*n) -- combined.M4 += 6.0*delta2 * (a.n*a.n*b.M2 + b.n*b.n*a.M2)/(combined.n*combined.n) + - + 6.0*delta2 * (and*and*bm2 + bnd*bnd*am2) / fi n2 + + 6.0*delta2 * (nda*nda*bm2 + ndb*ndb*am2) / fi n2 -- 4.0*delta*(a.n*b.M3 - b.n*a.M3) / combined.n; - + 4.0*delta*(and*bm3 - bnd*am3) / nd + + 4.0*delta*(nda*bm3 - ndb*am3) / nd -- | Efficiently compute the -- __length, mean, variance, skewness and kurtosis__ with a single pass. @@ -436,6 +450,7 @@ fastLMVSKu = getLMVSKu <$> foldLMVSKState {-# INLINE lmvsk0 #-} +lmvsk0 :: LMVSK lmvsk0 = LMVSK 0 0 0 0 0 -- | Performs the heavy lifting of fastLMVSK. This is exposed @@ -513,6 +528,8 @@ ,lrrYStats :: {-# UNPACK #-}!LMVSK } deriving (Show, Eq) +-- | The number of elements which make up this 'LinRegResult' +-- /Since: 0.1.4.1/ lrrCount :: LinRegResult -> Int lrrCount = lmvskCount . lrrXStats @@ -544,8 +561,8 @@ -} instance Semigroup LinRegState where {-# INLINE (<>) #-} - (LinRegState ax@(LMVSKState ax') ay@(LMVSKState ay') a_xy) - <> (LinRegState bx@(LMVSKState bx') by@(LMVSKState by') b_xy) + (LinRegState ax@(LMVSKState ax') ay a_xy) + <> (LinRegState bx@(LMVSKState bx') by b_xy) = LinRegState x y s_xy where an = lmvskCount ax' bn = lmvskCount bx' @@ -616,7 +633,7 @@ -- /Since: 0.1.4.0/ {-# INLINE getLinRegResult #-} getLinRegResult :: LinRegState -> LinRegResult -getLinRegResult (LinRegState vx@(LMVSKState vx') vy@(LMVSKState vy') s_xy) = LinRegResult slope intercept correlation statsx statsy where +getLinRegResult (LinRegState vx@(LMVSKState vx') vy s_xy) = LinRegResult slope intercept correl statsx statsy where n = lmvskCount vx' ndm1 = fromIntegral (n-1) -- slope = S_xy / (x_stats.Variance()*(n - 1.0)); @@ -625,7 +642,7 @@ slope = s_xy / lmvskVariance vx' intercept = yMean - slope*xMean t = sqrt xVar * sqrt yVar -- stddev x * stddev y - correlation = s_xy / (ndm1 * t) + correl = s_xy / (ndm1 * t) -- Need unbiased variance or correlation may be > ±1 statsx@(LMVSK _ xMean xVar _ _) = getLMVSKu vx statsy@(LMVSK _ yMean yVar _ _) = getLMVSKu vy @@ -639,7 +656,7 @@ {-# INLINE foldLinRegState #-} foldLinRegState :: Fold (Double,Double) LinRegState foldLinRegState = Fold step (LinRegState (LMVSKState lmvsk0) (LMVSKState lmvsk0) 0) id where - step st@(LinRegState vx@(LMVSKState vx') vy@(LMVSKState vy') s_xy) (x,y) = LinRegState vx2 vy2 s_xy' where + step (LinRegState vx@(LMVSKState vx') vy s_xy) (x,y) = LinRegState vx2 vy2 s_xy' where n = lmvskCount vx' nd = fromIntegral n nd1 = fromIntegral (n+1) @@ -656,7 +673,7 @@ -- returned by `fastLinearReg` 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) + 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) diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,25 +1,35 @@ +{-# LANGUAGE CPP #-} -import Test.Tasty +import Test.Tasty -- import Test.Tasty.SmallCheck as SC -import qualified Test.Tasty.QuickCheck as QC -import Test.Tasty.QuickCheck ((==>)) +import Test.Tasty.QuickCheck ((==>)) +import qualified Test.Tasty.QuickCheck as QC -import qualified Control.Foldl as F -import Control.Foldl.Statistics hiding (length) +#if MIN_VERSION_foldl(1,2,2) +import qualified Control.Foldl as F hiding (mean, variance) +#else +import qualified Control.Foldl as F +#endif -import qualified Data.Vector.Unboxed as U -import Test.QuickCheck.Instances +import Control.Foldl.Statistics hiding (length) -import qualified Statistics.Sample as S -import Statistics.Function (within) +import qualified Data.Vector.Unboxed as U +import Test.QuickCheck.Instances () -import Data.Profunctor +import qualified Statistics.Sample as S -import Data.Function (on) +import Data.Profunctor -import Data.Semigroup ((<>)) +import Data.Function (on) +import Data.Semigroup ((<>)) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +import Data.Monoid (mappend) +#endif + + toV :: [Double] -> U.Vector Double toV = U.fromList @@ -39,9 +49,10 @@ <*> skewness m <*> kurtosis m - +precision :: Double precision = 0.0000000001 +cmpLMVSK :: Double -> LMVSK -> LMVSK -> Bool cmpLMVSK prec a b = let t f = on (withinPCT prec) f a b in t lmvskMean @@ -76,6 +87,8 @@ , onVec "fastStdDev" $ \vec -> not (U.null vec) ==> F.fold fastStdDev (U.toList vec) == S.fastStdDev vec , let + -- TODO: Known failure when using + -- --quickcheck-replay '39 TFGenR A6EB566E901D554AAA13826C088B8831192E813D893D082A85F8A27C86D569E0 0 65535 16 0' in onVec ("fastLMVSK within " ++ show precision ++ " %") $ \vec -> U.length vec > 2 ==> let m = F.fold mean $ U.toList vec @@ -178,4 +191,4 @@ withinPCT :: Double -> Double -> Double -> Bool -withinPCT pct a b = abs (a-b) * 100 / (min `on` abs) a b < pct +withinPCT pct a b = abs (a - b) * 100 / (min `on` abs) a b < pct