foldl-statistics 0.1.0.0 → 0.1.1.0

raw patch · 3 files changed

+228/−104 lines, 3 filesdep ~math-functionsPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: math-functions

API changes (from Hackage documentation)

Files

• bench/Main.hs +109/−100
• foldl-statistics.cabal +2/−2
• src/Control/Foldl/Statistics.hs +117/−2

@@ -15,6 +15,9 @@ sample :: U.Vector Double sample = runST $ flip uniformVector 10000 =<< create +sample2 :: U.Vector (Double,Double)+sample2 = runST $ flip uniformVector 10000 =<< create+ absSample = U.map abs sample  -- Weighted test sample@@ -30,106 +33,112 @@  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 "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 "fastLMVSK"+                       -- T4 is strict in all arguments, so WHNF ok here+        [bench "C.F.Statistics"      $ whnf (\vec -> F.fold fastLMVSK (U.toList vec)) sample+        ]+      , bgroup "fastLinearReg"+        [bench "fastLinearReg"       $ whnf (\vec -> F.fold fastLinearReg (U.toList vec)) sample2+        ]+      ] -        , 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-                ]-            ]+    , 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+        ]+      ]+    ]

@@ -1,5 +1,5 @@ name:                foldl-statistics-version:             0.1.0.0+version:             0.1.1.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@@ -26,7 +26,7 @@   default-language:    Haskell2010   build-depends:       base >= 4.7 && < 5                        , foldl >= 1.1 && < 1.3-                       , math-functions >= 0.1 && < 0.2+                       , math-functions >= 0.1 && < 0.3                        , profunctors >= 5.2 && < 5.3  test-suite foldl-statistics-test

@@ -44,8 +44,12 @@     , fastVariance     , fastVarianceUnbiased     , fastStdDev+    , fastLMVSK+    , Stats4(..)+    , fastLinearReg+    , LinRegResult(..) -    -- $correlation+     , correlation      -- * References@@ -349,9 +353,117 @@ fastStdDev = sqrt fastVariance  --- $correlation++-- | When returned by `fastLMVSK`, contains the count, mean,+--  variance, skewness and kurtosis of a series of samples. --+-- _Since: 0.1.1.0_+data Stats4  = Stats4+  { stats4Count    :: {-# UNPACK #-}!Int+  , stats4Mean     :: {-# UNPACK #-}!Double+  , stats4Variance :: {-# UNPACK #-}!Double+  , stats4Skewness :: {-# UNPACK #-}!Double+  , stats4Kurtosis :: {-# UNPACK #-}!Double+  } deriving (Show, Eq)++-- | Efficiently compute the+-- __length, mean, variance, skewness and kurtosis__ with a single pass. --+-- _Since: 0.1.1.0_+{-# INLINE fastLMVSK #-}+fastLMVSK :: Fold Double Stats4+fastLMVSK = finalStats4 <$> foldStats4+++{-# INLINE stats40 #-}+stats40 = Stats4 0 0 0 0 0++-- This performs the grunt work of the fastLMVSK function above.+-- Note: The Stats4 returned by this doesn't contain the actual statistics+-- you're likely after, you must apply `finalStats4' to compute those.+--+-- See details on John Cook's article in the references below for+-- details.+{-# INLINE foldStats4 #-}+foldStats4 :: Fold Double Stats4+foldStats4 = Fold stepStats4 stats40 id++{-# INLINE stepStats4 #-}+stepStats4 :: Stats4 -> Double -> Stats4+stepStats4 (Stats4 n1 m1 m2 m3 m4) x = Stats4 n' m1' m2' m3' m4' where+  n' = n1+1+  delta = x - m1+  delta_n = delta / fromIntegral n'+  delta_n2 = delta_n * delta_n+  term1 = delta * delta_n * fromIntegral n1+  m1' = m1 + delta_n+  m4' = m4 + term1 * delta_n2 * fromIntegral (n'*n' - 3*n' + 3) + 6 * delta_n2 * m2 - 4 * delta_n * m3+  m3' = m3 + term1 * delta_n  * fromIntegral (n' - 2)           - 3 * delta_n  * m2+  m2' = m2 + term1+finalStats4 :: Stats4 -> Stats4+finalStats4 (Stats4 n m1 m2 m3 m4) = Stats4 n m1 m2' m3' m4' where+  nd = fromIntegral n+  m2' = m2 / (nd-1)+  m3' = sqrt nd * m3 * (m2 ** (-1.5))+  m4' = nd*m4 / (m2*m2) - 3.0++-- | When returned by `fastLinearReg`, contains the count,+--   slope, intercept and correlation of combining @(x,y)@ pairs.+--+-- _Since: 0.1.1.0_+data LinRegResult = LinRegResult+  {lrrCount       :: {-# UNPACK #-}!Int+  ,lrrSlope       :: {-# UNPACK #-}!Double+  ,lrrIntercept   :: {-# UNPACK #-}!Double+  ,lrrCorrelation :: {-# UNPACK #-}!Double+  } deriving (Show, Eq)++-- | Computes the __count, slope, (Y) intercept and correlation__ of @(x,y)@+--   pairs.+--+-- >>> F.fold fastLinearReg $ map (\x -> (x,3*x+7)) [1..100]+-- LinRegResult {lrrCount = 100, lrrSlope = 3.0,+--               lrrIntercept = 7.0, lrrCorrelation = 1.0}+--+-- >>> F.fold fastLinearReg $ map (\x -> (x,0.005*x*x+3*x+7)) [1..100]+-- LinRegResult {+--    lrrCount = 100,+--    lrrSlope = 3.5049999999999994,+--    lrrIntercept = -1.5849999999999795,+--    lrrCorrelation = 0.9993226275740273}+--+-- _Since: 0.1.1.0_+{-# INLINE fastLinearReg #-}+fastLinearReg :: Fold (Double,Double) LinRegResult+fastLinearReg = Fold step (V2 0 (V 0 0) (V 0 0) 0) final where+  step (V2 n v1@(V xMean xVar) v2@(V yMean _) s_xy) (x,y) = V2 (n+1) v1' v2' s_xy' where+    nd = fromIntegral n+    nd1 = fromIntegral (n+1)+    s_xy' = s_xy + (xMean - x)*(yMean - y)*nd/nd1+    v1' = stepV v1 n x+    v2' = stepV v2 n y+  final (V2 n v1@(V xMean xVar) v2@(V yMean yVar)  s_xy) = LinRegResult n slope intercept correlation where+    ndm1 = fromIntegral (n-1)+    slope = s_xy / xVar+    intercept = yMean - slope*xMean+    t = sqrt (xVar/ndm1) * sqrt (yVar/ndm1); -- stddev x * stddev y+    correlation = s_xy / (ndm1 * t)++data V2 = V2 {-# UNPACK #-}!Int {-# UNPACK #-}!V {-# UNPACK #-}!V {-# UNPACK #-}!Double++{-# INLINE stepV #-}+stepV :: V -> Int -> Double -> V+stepV (V m1 m2) n1 x = V m1' m2' where+  delta = x - m1+  delta_n = delta / fromIntegral (n1+1)+  term1 = delta * delta_n * fromIntegral n1+  m1' = m1 + delta_n+  m2' = m2 + term1++++-- | Given the mean and standard deviation of two distributions, computes+--   the correlation between them. 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)@@ -376,6 +488,9 @@ -- * West, D.H.D. (1979) Updating mean and variance estimates: an --   improved method. /Communications of the ACM/ --   22(9):532&#8211;535. <http://doi.acm.org/10.1145/359146.359153>+--+-- * John D. Cook. Computing skewness and kurtosis in one pass+--   <http://www.johndcook.com/blog/skewness_kurtosis/>