diff --git a/bench/Main.hs b/bench/Main.hs --- a/bench/Main.hs +++ b/bench/Main.hs @@ -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 + ] + ] + ] diff --git a/foldl-statistics.cabal b/foldl-statistics.cabal --- a/foldl-statistics.cabal +++ b/foldl-statistics.cabal @@ -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 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 @@ -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–535. +-- +-- * John D. Cook. Computing skewness and kurtosis in one pass +--