Provide approximate equivalents for non fusing functions

[mchav]

A number of functions that rely on robust sum don't fuse because of the constraint that the sum truncate as little as possible (using Numeric.Sum). I haven't tested out how much accuracy is lost but I think it might be worth having a faster version of functions like mean that does fuse by folding over the elements, accumulating the sum and length in one pass then finally adding.

Maybe called approxMean?

[mchav]

Similarly it would be good to have generic versions that operate on different floating point representations e.g float16 and float32 instead of being locked in to using the highest precision float representation.

[Shimuuar]

Kahan summation should cause any problem. It's plain fold and stream fusion has no problem with it. Using length breaks fusion since it requires materialized vector.

Float 32/16 is a bit tricky. Frist it could cause performance problems: making field polymorphic means it couldn't be unboxed and we're completely at mercy of optimizer. It may or may not unbox it. Also there' little reason to use smaller floating point numbers as an accumulators when working on CPU. They aren't faster (unless SIMD is used) and much more likely to run into numeric problems

P.S. you may use monoid-statisitcs which provides statistics like mean in form of pure folds.

[mchav]

I see. Maybe my diagnosis was premature then - I misread and thought we were using sumF.

Why doesn't mean fuse then - according to this comment?

-- | /O(n)/ Arithmetic mean.  This uses Kahan-Babuška-Neumaier
-- summation, so is more accurate than 'welfordMean' unless the input
-- values are very large. This function is not subject to stream
-- fusion.
mean :: (G.Vector v Double) => v Double -> Double
mean xs = sum xs / fromIntegral (G.length xs)
{-# SPECIALIZE mean :: U.Vector Double -> Double #-}
{-# SPECIALIZE mean :: V.Vector Double -> Double #-}

[mchav]

Oh. I just read the comment above that says:

-- Because of the need for two passes, these functions are /not/
-- subject to stream fusion.

Would it make sense to compute mean in one pass like so:

let (total, n) = G.foldl' (\(total, n) v -> (total + v, n +1) (0, 0) samp
in total / n

[Shimuuar]

It's just very poorly worded. Very poorly. In fact there're two reasons:

1. xs is used twice. That's enough to break fusion
2. length is not fusible.

Yes that will work. Except tuple is lazy in both argument and this implementation may leak despite using foldl' Forcing tuple doesn't force tuple elements

[mchav]

Did some benchmarking to see the actual difference:

{-# LANGUAGE NumericUnderscores #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE BangPatterns #-}

import qualified Data.Vector.Unboxed as VU

import Control.Monad (replicateM)
import Criterion.Main
import System.Random (randomRIO)

import qualified Data.Vector.Generic as G
import qualified Statistics.Sample as S

custom :: VU.Vector Double -> IO ()
custom samp = let
    (!total, !n) = G.foldl' (\(!total, !n) v -> (total + v, n + 1))  (0 :: Double, 0 :: Int) samp
  in print $ total / fromIntegral n

outOfBox :: VU.Vector Double -> IO ()
outOfBox samp = print $ S.mean samp 

main = do
  ns <- VU.replicateM 300_000_000 (randomRIO (-20.0 :: Double, 20.0))
  defaultMain [
    bgroup "stats" [ bench     "custom" $ nfIO (custom ns)
                   , bench   "outOfBox" $ nfIO (outOfBox ns)
                   ]
    ]

The results are as follows:

Current implementation:

time                 1.389 s    (1.333 s .. 1.492 s)
                     0.999 R²   (0.998 R² .. 1.000 R²)
mean                 1.355 s    (1.343 s .. 1.374 s)
std dev              18.43 ms   (1.868 ms .. 23.20 ms)
variance introduced by outliers: 19% (moderately inflated)

Accumulator implementation:

time                 438.6 ms   (415.4 ms .. 481.9 ms)
                     0.999 R²   (0.997 R² .. 1.000 R²)
mean                 451.5 ms   (445.9 ms .. 461.3 ms)
std dev              9.627 ms   (2.839 ms .. 12.87 ms)
variance introduced by outliers: 19% (moderately inflated)

Seems to be a sizeable speedup.

[mchav]

Even when I try welfordMean (which is a single pass mean) I get a pretty slow impementation:

time                 1.932 s    (1.828 s .. 2.048 s)
                     1.000 R²   (0.998 R² .. 1.000 R²)
mean                 1.913 s    (1.896 s .. 1.936 s)
std dev              22.91 ms   (7.176 ms .. 30.24 ms)
variance introduced by outliers: 19% (moderately inflated)

Can you reproduce a similar benchmarking trend?

[axman6]

I believe a lot of these problems are solved by my foldl-statistics package, which implements all of Statistics.Sample (at least how it existed a few years ago) using the foldl package. The implementation of the efficient mean above would simply be:

import Control.Foldl as Fold

fastMean :: Vector Double -> Double
fastMean = Foldl.fold (sum / genericLength)

I'd love to see my package retired and many of the algorithms in Statistics.Sample implemented as Folds and exposed. The performance last time I benchmarked showed my packages was equal to or slightly faster than statistics back in the day, but you get significantly more flexibility by being able to combine Folds applicatively to make arbitrarily complex single-pass computations.

(@mchav I pinged you on the data-haskell chat but it looks pretty dead, so moving the convo here).

[mchav]

Have you profiled this approach recently? I know this library ends up piggy backing off of improvements in vector.

[axman6]

I haven't, but those improvements should be applicable to this approach too - the actual fold that would be used for most of the functions which currently use G.foldl would become (something like this, I'm unsure which inline and specialise pragmas would be best):

ffold :: G.Vector v a => F.Fold a b -> v a -> b
ffold f v = F.purely_ (\s i -> G.foldl s i v) f
{-# INLINE ffold #-}
{-# SPECIALIZE ffold :: F.Fold Double b -> U.Vector Double -> b #-}
{-# SPECIALIZE ffold :: F.Fold Double b -> V.Vector Double -> b #-}

I'll make an initial attempt and see how it goes.

[mchav]

Thanks @axman6 I'll hot swap it in a branch and see what performance looks like!

[axman6]

I ran some initial benchmarks, and the performance appears to be very similar (my machine is definitely not an ideal benchmark environment at the moment, plenty of other thinds going on in the background):

                           , master   , foldl    , % of master, master      , foldl
Name                       , Mean (ps), Mean (ps),            , 2*Stdev (ps), 2*Stdev (ps)
All.sample.range           , 8216934  , 8410186  ,      102.35, 516844      , 838116
All.sample.mean            , 28347814 , 26064050 ,       91.94, 1834328     , 2164224
All.sample.meanWeighted    , 52317871 , 50228173 ,       96.00, 4359130     , 4910178
All.sample.harmonicMean    , 6512393  , 6546374  ,      100.52, 430016      , 536538
All.sample.geometricMean   , 48041699 , 48252929 ,      100.43, 3829484     , 3506958
All.sample.variance        , 42851293 , 42593981 ,       99.39, 3917158     , 2085504
All.sample.varianceUnbiased, 42557666 , 42652026 ,      100.22, 3588736     , 3757096
All.sample.varianceWeighted, 60825317 , 60919604 ,      100.15, 4340024     , 5591370
All.sample.stdDev          , 42590405 , 42733349 ,      100.33, 3409676     , 3996206
All.sample.skewness        , 112154589, 113525439,      101.22, 7980622     , 9072452
All.sample.kurtosis        , 124605078, 119529199,       95.92, 8816308     , 9309160
All.sample.C.M. 2          , 71077685 , 73399584 ,      103.26, 5378790     , 4276260
All.sample.C.M. 3          , 79193408 , 79100683 ,       99.88, 7121570     , 7726006
All.sample.C.M. 4          , 85476806 , 86221044 ,      100.87, 6667806     , 7727098
All.sample.C.M. 5          , 96715087 , 95239404 ,       98.47, 8259694     , 6979196
-- Not implemented yet
All.sample.pearson         , 99972656 , 100695898,      100.72, 7165680     , 7674690
All.sample.covariance      , 71694995 , 67749365 ,       94.49, 5445392     , 6700342
All.sample.correlation     , 99785937 , 99124438 ,       99.33, 7342992     , 7913474
All.sample.covariance2     , 70288110 , 68865185 ,       97.97, 5051024     , 6739430
All.sample.correlation2    , 98770507 , 98927978 ,      100.15, 8162520     , 7703952

I'll push up my (very much WIP) branch for you.

[axman6]

https://github.com/haskell/statistics/compare/master...axman6:statistics:foldl-statistics?expand=1

[Shimuuar]

@mchav your implementation uses plain double summation while Statistics.Sample.mean uses compensated summation (effectively working in quad precision at cost of 2× computation). Switching fold to KBN summation gives same performance