foldl-statistics 0.1.3.0 → 0.1.4.0

raw patch · 5 files changed

+171/−59 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.md +3/−0
• bench/Main.hs +4/−4
• foldl-statistics.cabal +1/−1
• src/Control/Foldl/Statistics.hs +135/−36
• test/Spec.hs +28/−18

@@ -1,3 +1,6 @@+# 0.1.4.0+- Added monoidal interface to linear regression+ # 0.1.3.0 - Added unbiased versions of LMVSK functions 

@@ -74,11 +74,11 @@         ]       , 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-        ]+          [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-        ]+          [bench "fastLinearReg"       $ whnf (\vec -> F.fold fastLinearReg (U.toList vec)) sample2+          ]       ]      , bgroup "requiring the mean"

@@ -1,5 +1,5 @@ name:                foldl-statistics-version:             0.1.3.0+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

@@ -51,10 +51,13 @@     , foldLMVSKState     , getLMVSK     , getLMVSKu++    -- ** Linear Regression     , fastLinearReg+    , foldLinRegState+    , getLinRegResult     , LinRegResult(..)--+    , LinRegState     , correlation      -- * References@@ -503,58 +506,154 @@ -- -- /Since: 0.1.1.0/ data LinRegResult = LinRegResult-  {lrrCount       :: {-# UNPACK #-}!Int-  ,lrrSlope       :: {-# UNPACK #-}!Double+  {lrrSlope       :: {-# UNPACK #-}!Double   ,lrrIntercept   :: {-# UNPACK #-}!Double   ,lrrCorrelation :: {-# UNPACK #-}!Double+  ,lrrXStats      :: {-# UNPACK #-}!LMVSK+  ,lrrYStats      :: {-# UNPACK #-}!LMVSK   } deriving (Show, Eq) --- | Computes the __count, slope, (Y) intercept and correlation__ of @(x,y)@---   pairs.+lrrCount :: LinRegResult -> Int+lrrCount = lmvskCount . lrrXStats++-- | The Monoidal state used to compute linear regression, see `fastLinearReg`. --+-- /Since: 0.1.4.0/+data LinRegState = LinRegState+  {-# UNPACK #-}!LMVSKState+  {-# UNPACK #-}!LMVSKState+  {-# UNPACK #-}!Double+++{-+RunningRegression operator+(const RunningRegression a, const RunningRegression b)+{+    RunningRegression combined;++    combined.x_stats = a.x_stats + b.x_stats;+    combined.y_stats = a.y_stats + b.y_stats;+    combined.n = a.n + b.n;++    double delta_x = b.x_stats.Mean() - a.x_stats.Mean();+    double delta_y = b.y_stats.Mean() - a.y_stats.Mean();+    combined.S_xy = a.S_xy + b.S_xy ++    double(a.n*b.n)*delta_x*delta_y/double(combined.n);++    return combined;+}+-}+instance Semigroup LinRegState where+  {-# INLINE (<>) #-}+  (LinRegState ax@(LMVSKState ax') ay@(LMVSKState ay') a_xy)+   <> (LinRegState bx@(LMVSKState bx') by@(LMVSKState by') b_xy)+   = LinRegState x y s_xy where+    an = lmvskCount ax'+    bn = lmvskCount bx'+    x = ax <> bx+    y = ay <> by+    delta_x = lmvskMean (getLMVSK bx) - lmvskMean (getLMVSK ax)+    delta_y = lmvskMean (getLMVSK by) - lmvskMean (getLMVSK ay)+    s_xy = a_xy+b_xy + fromIntegral (an*bn) * delta_x * delta_y/fromIntegral (an+bn)+++instance Monoid LinRegState where+  {-# INLINE mempty #-}+  mempty = LinRegState mempty mempty 0+  {-# INLINE mappend #-}+  mappend = (<>)++++-- | Computes the __slope, (Y) intercept and correlation__ of @(x,y)@+--   pairs, as well as the `LMVSK` stats for both the x and y series.+-- -- >>> F.fold fastLinearReg $ map (\x -> (x,3*x+7)) [1..100]--- LinRegResult {lrrCount = 100, lrrSlope = 3.0,---               lrrIntercept = 7.0, lrrCorrelation = 1.0}+-- LinRegResult+--   {lrrSlope = 3.0+--   , lrrIntercept = 7.0+--   , lrrCorrelation = 100.0+--   , lrrXStats = LMVSK+--       {lmvskCount = 100+--       , lmvskMean = 50.5+--       , lmvskVariance = 833.25+--       , lmvskSkewness = 0.0+--       , lmvskKurtosis = -1.2002400240024003}+--   , lrrYStats = LMVSK+--       {lmvskCount = 100+--       , lmvskMean = 158.5+--       , lmvskVariance = 7499.25+--       , lmvskSkewness = 0.0+--       , lmvskKurtosis = -1.2002400240024003}+--   } -- -- >>> 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}+-- LinRegResult+--   {lrrSlope = 3.5049999999999994+--   , lrrIntercept = -1.5849999999999795+--   , lrrCorrelation = 99.93226275740273+--   , lrrXStats = LMVSK+--       {lmvskCount = 100+--       , lmvskMean = 50.5+--       , lmvskVariance = 833.25+--       , lmvskSkewness = 0.0+--       , lmvskKurtosis = -1.2002400240024003}+--   , lrrYStats = LMVSK+--       {lmvskCount = 100+--       , lmvskMean = 175.4175+--       , lmvskVariance = 10250.37902625+--       , lmvskSkewness = 9.862971188165422e-2+--       , lmvskKurtosis = -1.1923628437011482}+--   } -- -- /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)+fastLinearReg = getLinRegResult <$> foldLinRegState -data V2 = V2 {-# UNPACK #-}!Int {-# UNPACK #-}!V {-# UNPACK #-}!V {-# UNPACK #-}!Double+-- | Produces the slope, Y intercept, correlation and LMVSK stats from a+--   `LinRegState`.+--+-- /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+  n                               = lmvskCount vx'+  ndm1                            = fromIntegral (n-1)+  -- slope = S_xy / (x_stats.Variance()*(n - 1.0));+  -- in LMVSKState, 'lmvskVariance' hasn't been divided+  -- by (n-1), so division not necessary+  slope                           = s_xy / lmvskVariance vx'+  intercept                       = yMean - slope*xMean+  t                               = sqrt xVar * sqrt yVar -- stddev x * stddev y+  correlation                     = s_xy / (ndm1 * t)+  -- Need unbiased variance or correlation may be > ±1+  statsx@(LMVSK _ xMean xVar _ _) = getLMVSKu vx+  statsy@(LMVSK _ yMean yVar _ _) = getLMVSKu vy -{-# 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 +-- | Performs the heavy lifting for `fastLinReg`. Exposed because `LinRegState`+--  is a `Monoid`, allowing statistics to be computed on datasets in parallel+--  and combined afterwards.+--+-- /Since: 0.1.4.0/+{-# 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+    n     = lmvskCount vx'+    nd    = fromIntegral n+    nd1   = fromIntegral (n+1)+    s_xy' = s_xy + (xMean - x)*(yMean - y)*nd/nd1+    xMean = lmvskMean (getLMVSK vx)+    yMean = lmvskMean (getLMVSK vy)+    vx2   = stepLMVSKState vx x+    vy2   = stepLMVSKState vy y   -- | Given the mean and standard deviation of two distributions, computes---   the correlation between them.+--   the correlation between them, given the means and standard deviation+--   of the @x@ and @y@ series. The results may be more accurate than those+--   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)

@@ -39,6 +39,17 @@   <*> skewness m   <*> kurtosis m ++precision = 0.0000000001++cmpLMVSK prec a b = let+  t f = on (withinPCT prec) f a b+  in t lmvskMean+     && t lmvskVariance+     && t lmvskKurtosis+     && t lmvskSkewness+     && ((==) `on` lmvskCount) a b+ main :: IO () main = defaultMain $     testGroup "Results match Statistics.Sample"@@ -58,15 +69,7 @@                 ]              , testGroup "Single-pass functions" $-                let precision = 0.0000000001-                    cmp prec a b = let-                      t f = on (withinPCT prec) f a b-                      in t lmvskMean-                         && t lmvskVariance-                         && t lmvskKurtosis-                         && t lmvskSkewness-                         && ((==) `on` lmvskCount) a b-                in [ onVec "fastVariance" $ \vec ->+                [ 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@@ -78,12 +81,12 @@                       m         = F.fold mean $ U.toList vec                       fast      = F.fold fastLMVSK $ U.toList vec                       reference = F.fold (testLMVSK m) $ U.toList vec-                      in cmp precision fast reference+                      in cmpLMVSK precision fast reference                 , QC.testProperty "LMVSKSemigroup" $ \v1 v2 ->                     U.length v1 > 2 && U.length v2 > 2 && U.sum (mappend v1 v1) /= U.product (mappend v1 v1) ==> let                       sep = getLMVSK $ F.fold foldLMVSKState (U.toList v1) <> F.fold foldLMVSKState (U.toList v2)                       tog = F.fold fastLMVSK (U.toList v1 ++ U.toList v2)-                      in cmp precision sep tog+                      in cmpLMVSK precision sep tog                         || isNaN (lmvskKurtosis sep)                         || isNaN (lmvskKurtosis tog)                 ]@@ -146,19 +149,26 @@                         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)+                    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))+                        (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-+                , QC.testProperty "LinRegState Semigroup" $ \v1 v2 ->+                    U.length v1 > 2 && U.length v2 > 2+                    && U.sum (U.map fst (mappend v1 v1)) /= U.product (U.map fst (mappend v1 v1))+                    && U.sum (U.map snd (mappend v1 v1)) /= U.product (U.map snd (mappend v1 v1)) ==> let+                      sep = getLinRegResult $ F.fold foldLinRegState (U.toList v1) <> F.fold foldLinRegState (U.toList v2)+                      tog = F.fold fastLinearReg (U.toList v1 ++ U.toList v2)+                      in (cmpLMVSK precision (lrrXStats sep) (lrrXStats tog)+                         && cmpLMVSK precision (lrrYStats sep) (lrrYStats tog))+                        || isNaN (lmvskKurtosis (lrrXStats sep))+                        || isNaN (lmvskKurtosis (lrrYStats sep))+                        || isNaN (lmvskKurtosis (lrrXStats tog))+                        || isNaN (lmvskKurtosis (lrrYStats tog))                 ]             ]         ]