diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,6 @@ +# 0.1.4.0 +- Added monoidal interface to linear regression + # 0.1.3.0 - Added unbiased versions of LMVSK functions diff --git a/bench/Main.hs b/bench/Main.hs --- a/bench/Main.hs +++ b/bench/Main.hs @@ -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" 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.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 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 @@ -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) diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -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)) ] ] ]