Folding Incremental Averages in Haskell

In a recent post, SFTank describes an interesting application of foldr to computing the mean of a list of numbers via a fold. Basically, their approach -- rather than summing up -- computes the mean m in an incremental fashion via

m_n+1 = m_n + (x_n+1 - m_n) / (n + 1) (1)

which we can write in Haskell as


> increment :: (Fractional a, Integral b) => b -> a -> a -> a
> increment n mean x = mean + (x - mean) / fromIntegral n

The post then proceeds to define an ``index-aware'' foldr (i.e. a foldr where the folding knows the list index of the element it's folding over). Let's factor that out:


> foldri :: Integral n => (a -> b -> n -> b) -> b -> [a] -> b
> foldri f y_0 l = fst $ foldr fi (y_0, 1) l
>     where 
>         fi x (y, i) = (f x y i, i + 1)

I think I can see some monadic structure shining through, but we'll ignore that for now. Instead, let's first proceed to implement their averaging function


> average l = foldri inc 0
>     where
>         inc x m i = increment i m x

This is only semi-elegant, since we're basically just unnecessarily shuffling around arguments, but bear with me, since it's not what I'm after anyhow. Instead, consider this: the increment not only works for obtaining the average of x, (i'll call it <x>), but (obviously) also the average of the square of x, and in fact any <xⁿ> (i.e. moments about the origin). Why would I think this is cool? Well, look at the form of e.g. the sample variance:

σ_x = <x²> - <x>² (2)

so in principle, we can obtain the sample mean and variances in one pass over the data. In fact, in a lazy language, we can obtain all <xⁿ> with just one pass, and then conveniently compute properties such as mean, variance, skew and kurtosis. Which is exactly what I'll do.


> powers :: Num a => a -> [a]
> powers x = map (\n -> x^n) [1..]

> moments' :: (Fractional a, Integral i) => a -> [a] -> i -> [a]
> moments' x m n = zipWith (increment n) m (powers x)

> moments = foldri moments' (repeat 0)

So now, we get a whole bunch of stuff gratis

*Main> take 5 $ moments [1, 2, 3]
[2.0,4.666666666666667,12.0,32.666666666666664,92.0]

I like this. Next we can simply define the mean as


> mean = head . moments

Here's some second-order stuff,


> variance x = let [mu, sigma] = take 2 $ moments x
>    in
>        sigma - mu^2

> stddev = sqrt . variance

For skew and kurtosis, we'll need moments about the mean, rather than the moments about the origin. While that would be easy to implement in a general fashion (and a good excuse for doing factorials and binomial coefficients in haskell -- in spite of what your friends might say) I will not do that here, and instead just hand-code the definitions from wikipedia. Being hand-coded and not properly tested, they're probably buggy, but I assume by now you get the gist:


> skew x = let 
>     [mu_1', mu_2', mu_3'] = take 3 $ moments x
>     mu_3 = mu_3' - 3 * mu_1' * mu_2' + 2 * (mu_1' ^ 3)
>     sigma = sqrt (mu_2' - (mu_1' ^ 2))
>     in
>         mu_3 / (sigma ^ 3)
>
> kurtosis x = let
>     [mu_1', mu_2', mu_3', mu_4'] = take 4 $ moments x
>     sigma2 = mu_2' - mu_1'
>     mu_4 = mu_4' - 4 * mu_1' * mu_3' + 6 * (mu_1'^2) * mu_2' - 3 * (mu_1'^4)
>     in
>       mu_4 / (sigma2^2)