diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,7 @@ +# 0.1.1.0 +- Add fastLMVSK (length, mean, variance, skewness and kurtosis) +- Add fastLinearReg (count, slope, (Y) intercept and correlation of `(x,y)`) + + +# 0.1.0.0 +- Initial release diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,45 @@ +# foldl-statistics +A reimplementation of the [Statistics.Sample](https://hackage.haskell.org/package/statistics/docs/Statistics-Sample.html) +module using the [foldl](https://www.stackage.org/lts-5.1/package/foldl) package. +The intention of this package is to allow these algorithms to be used on a much broader set of data input types, +including lists and streaming libraries such as `conduit` and `pipes`, and any other type which is `Foldable`. + +All statistics in this package can be computed with no more than two passes over the data - once to compute the mean and once to compute +any statistics which require the mean. this is achieved because foldl `Fold`s are `Applicative`, which means that to compute, for example, the first 4 central moments as well as the count, the following could be used: + +```haskell +import Control.Foldl as F + +... + +dataseries :: [Double] +dataseries = ... + + ... + let m = F.fold mean dataseries + (c2,c3,c4,c5,n) = flip F.fold dataseries $ + (\(c2,c3) (c4,c5) n -> (c2,c3,c4,c5,n)) + <$> centralMoment 2 3 m + <*> centralMoment 4 5 m + <*> F.length +``` + +which traverses the data twice, once to compute the mean `m`, and once to compute all the central moments and the count concurrently. This brings along with it for free the ability to compute streaming statistics, such as the mean of all data seen so far, using the `foldl`'s `scan` function. + +Where possible, care has been taken to ensure the numerical stability of the computation of statistics. + +Several algorithms require the mean of the data to be known before computing the statistic, such as `skewness`, `kurtosis` and other `centralMoment`s. +There are 'fast' implementations for calculating the variance, unbiased variance and standard deviation, which can be computed without knowing the mean +*a priori*, but which may produce less accurate results. + +## Performance & Correctness +Benchmarks are included comparing performance to the [statistics](https://hackage.haskell.org/package/statistics) package. In nearly all cases, the implementations in this package perform better than those in `statistics` on the same inputs, and in several cases, performing two passes (to compute the mean and another statistic) is faster than the equivalent `statistics` implementation. + +This speed has not come at the cost of correctness; all `Fold`s are tested against their `statistics` counterparts to ensure the results are identical. + +These results can be confirmed by running + + stack build --test --bench --benchmark-arguments "--output bench.html" + +which will print out the results of the tests against `statistics` and then run the benchmark (this may take several minutes and is best run on a "quiet" machine which is doing very little other than running the benchmark). The results of the benchmarking are then available in the file `bench.html`. + diff --git a/bench/Main.hs b/bench/Main.hs --- a/bench/Main.hs +++ b/bench/Main.hs @@ -12,15 +12,19 @@ -- Test sample +{-# NOINLINE sample #-} sample :: U.Vector Double sample = runST $ flip uniformVector 10000 =<< create +{-# NOINLINE sample2 #-} sample2 :: U.Vector (Double,Double) sample2 = runST $ flip uniformVector 10000 =<< create +{-# NOINLINE absSample #-} absSample = U.map abs sample -- Weighted test sample +{-# NOINLINE sampleW #-} sampleW :: U.Vector (Double,Double) sampleW = U.zip sample (U.reverse sample) @@ -55,7 +59,7 @@ , bench "Statistics.Sample" $ nf S.geometricMean absSample ] ] - , bgroup "Single-pass functions" + , bgroup "Single-pass" [ bgroup "fastVariance" [ bench "C.F.Statistics" $ nf (\vec -> F.fold fastVariance (U.toList vec)) sample , bench "Statistics.Sample" $ nf S.fastVariance sample @@ -77,7 +81,7 @@ ] ] - , bgroup "Functions requiring the mean" + , bgroup "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 @@ -100,7 +104,7 @@ ] ] - , bgroup "Functions over central moments" + , bgroup "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 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.1.0 +version: 0.1.2.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 @@ -17,7 +17,7 @@ copyright: 2016 Data61 (CSIRO) category: Math, Statistics build-type: Simple --- extra-source-files: +extra-source-files: CHANGELOG.md, README.md cabal-version: >=1.10 library @@ -28,6 +28,7 @@ , foldl >= 1.1 && < 1.3 , math-functions >= 0.1 && < 0.3 , profunctors >= 5.2 && < 5.3 + , semigroups test-suite foldl-statistics-test type: exitcode-stdio-1.0 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 @@ -45,7 +45,10 @@ , fastVarianceUnbiased , fastStdDev , fastLMVSK - , Stats4(..) + , LMVSK(..) + , LMVSKState + , foldLMVSKState + , getLMVSK , fastLinearReg , LinRegResult(..) @@ -61,6 +64,7 @@ import Control.Foldl as F import qualified Control.Foldl import Data.Profunctor +import Data.Semigroup import Numeric.Sum (KBNSum, kbn, add, zero) @@ -357,60 +361,125 @@ -- | 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 +-- /Since: 0.1.1.0/ +data LMVSK = LMVSK + { lmvskCount :: {-# UNPACK #-}!Int + , lmvskMean :: {-# UNPACK #-}!Double + , lmvskVariance :: {-# UNPACK #-}!Double + , lmvskSkewness :: {-# UNPACK #-}!Double + , lmvskKurtosis :: {-# UNPACK #-}!Double } deriving (Show, Eq) +newtype LMVSKState = LMVSKState LMVSK + +instance Monoid LMVSKState where + {-# INLINE mempty #-} + mempty = LMVSKState lmvsk0 + {-# INLINE mappend #-} + mappend = (<>) + +instance Semigroup LMVSKState where + {-# INLINE (<>) #-} + (LMVSKState (LMVSK an am1 am2 am3 am4)) <> (LMVSKState (LMVSK bn bm1 bm2 bm3 bm4)) + = LMVSKState (LMVSK n m1 m2 m3 m4) where + fi :: Int -> Double + fi = fromIntegral + -- combined.n = a.n + b.n; + n = an+bn + n2 = n*n + nd = fi n + and = fi an + bnd = fi bn + -- delta = b.M1 - a.M1; + delta = bm1 - am1 + -- delta2 = delta*delta; + delta2 = delta*delta + -- delta3 = delta*delta2; + delta3 = delta*delta2 + -- delta4 = delta2*delta2; + delta4 = delta2*delta2 + -- combined.M1 = (a.n*a.M1 + b.n*b.M1) / combined.n; + m1 = (and*am1 + bnd*bm1 ) / nd + -- combined.M2 = a.M2 + b.M2 + delta2*a.n*b.n / combined.n; + m2 = am2 + bm2 + delta2*and*bnd / 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 + -- combined.M3 += 3.0*delta * (a.n*b.M2 - b.n*a.M2) / combined.n; + + 3.0*delta * (and*bm2 - bnd*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) + -- 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 + -- 4.0*delta*(a.n*b.M3 - b.n*a.M3) / combined.n; + + 4.0*delta*(and*bm3 - bnd*am3) / nd + -- | Efficiently compute the -- __length, mean, variance, skewness and kurtosis__ with a single pass. -- --- _Since: 0.1.1.0_ +-- /Since: 0.1.1.0/ {-# INLINE fastLMVSK #-} -fastLMVSK :: Fold Double Stats4 -fastLMVSK = finalStats4 <$> foldStats4 +fastLMVSK :: Fold Double LMVSK +fastLMVSK = getLMVSK <$> foldLMVSKState -{-# INLINE stats40 #-} -stats40 = Stats4 0 0 0 0 0 +{-# INLINE lmvsk0 #-} +lmvsk0 = LMVSK 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. +-- | Performs the heavy lifting of fastLMVSK. This is exposed +-- because the internal `LMVSKState` is monoidal, allowing you +-- to run these statistics in parallel over datasets which are +-- split and then combine the results. -- --- See details on John Cook's article in the references below for --- details. -{-# INLINE foldStats4 #-} -foldStats4 :: Fold Double Stats4 -foldStats4 = Fold stepStats4 stats40 id +-- /Since: 0.1.2.0/ +{-# INLINE foldLMVSKState #-} +foldLMVSKState :: Fold Double LMVSKState +foldLMVSKState = Fold stepLMVSKState (LMVSKState lmvsk0) 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 +{-# INLINE stepLMVSKState #-} +stepLMVSKState :: LMVSKState -> Double -> LMVSKState +stepLMVSKState (LMVSKState (LMVSK n1 m1 m2 m3 m4)) x = LMVSKState $ LMVSK n m1' m2' m3' m4' where + fi :: Int -> Double + fi = fromIntegral + -- long long n1 = n; + -- n++; + n = n1+1 + -- delta = x - M1; + delta = x - m1 + -- delta_n = delta / n; + delta_n = delta / fi n + -- delta_n2 = delta_n * delta_n; + delta_n2 = delta_n * delta_n + -- term1 = delta * delta_n * n1; + term1 = delta * delta_n * fi n1 + -- M1 += delta_n; 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 + -- M4 += term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3; + m4' = m4 + term1 * delta_n2 * fi (n*n - 3*n + 3) + 6 * delta_n2 * m2 - 4 * delta_n * m3 + -- M3 += term1 * delta_n * (n - 2) - 3 * delta_n * M2; + m3' = m3 + term1 * delta_n * fi (n - 2) - 3 * delta_n * m2 + -- M2 += term1; m2' = m2 + term1 -finalStats4 :: Stats4 -> Stats4 -finalStats4 (Stats4 n m1 m2 m3 m4) = Stats4 n m1 m2' m3' m4' where + +-- | Returns the stats which have been computed in a LMVSKState. +-- +-- /Since: 0.1.2.0/ +getLMVSK :: LMVSKState -> LMVSK +getLMVSK (LMVSKState (LMVSK n m1 m2 m3 m4)) = LMVSK n m1 m2' m3' m4' where nd = fromIntegral n + -- M2/(n-1.0) m2' = m2 / (nd-1) - m3' = sqrt nd * m3 * (m2 ** (-1.5)) - m4' = nd*m4 / (m2*m2) - 3.0 + -- sqrt(double(n)) * M3/ pow(M2, 1.5) + m3' = sqrt nd * m3 / (m2 ** 1.5) + -- double(n)*M4 / (M2*M2) - 3.0 + 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_ +-- /Since: 0.1.1.0/ data LinRegResult = LinRegResult {lrrCount :: {-# UNPACK #-}!Int ,lrrSlope :: {-# UNPACK #-}!Double @@ -432,7 +501,7 @@ -- lrrIntercept = -1.5849999999999795, -- lrrCorrelation = 0.9993226275740273} -- --- _Since: 0.1.1.0_ +-- /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 diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -15,7 +15,11 @@ import Data.Profunctor +import Data.Function (on) +import Data.Semigroup ((<>)) + + toV :: [Double] -> U.Vector Double toV = U.fromList @@ -26,6 +30,15 @@ onVec2 :: String -> (U.Vector (Double,Double) -> QC.Property) -> TestTree onVec2 str f = QC.testProperty str (f . U.fromList) + +testLMVSK :: Double -> Fold Double LMVSK +testLMVSK m = LMVSK + <$> F.length + <*> mean + <*> varianceUnbiased m + <*> skewness m + <*> kurtosis m + main :: IO () main = defaultMain $ testGroup "Results match Statistics.Sample" @@ -44,13 +57,35 @@ in F.fold geometricMean (U.toList vec') == S.geometricMean vec' ] - , testGroup "Single-pass functions" - [ onVec "fastVariance" $ \vec -> + , 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 -> 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 , onVec "fastStdDev" $ \vec -> not (U.null vec) ==> F.fold fastStdDev (U.toList vec) == S.fastStdDev vec + , let + in onVec ("fastLMVSK within " ++ show precision ++ " %") $ \vec -> + U.length vec > 2 ==> let + 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 + , 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 + || isNaN (lmvskKurtosis sep) + || isNaN (lmvskKurtosis tog) ] ] @@ -130,3 +165,7 @@ between :: (Double,Double) -> Double -> Bool between (lo,hi) = \x -> lo <= x && x <= hi + + +withinPCT :: Double -> Double -> Double -> Bool +withinPCT pct a b = abs (a-b) * 100 / (min `on` abs) a b < pct