compute kth central moments in one pass

Continuing on the same topic as my previous post, its nice to be able to gather all the kth order moments in a single pass.

Last time I mentioned the boost/accumulators example, but you will have noticed two issues if you use that. Firstly, moment<k> tag will give you the kth simple moment relative to zero, whereas we often want the kth central moment of a sequence relative to the mean. Secondly, although boosts accumulator is well written it does seem to take a while to compile [~ 12 seconds for code using moment<12>].

After some playing around Ive got a faster simpler approach, where the inner loop accumulates kth powers of the element. After you’ve run the sequence through, you can then easily extract variance, and all the kth central moments. So in adding the more general case of kth moments, Ive made the particular variance case simpler. That often seems to happen in programming and math!

algebra

First a bit of math and then the code. We want to express the kth central moment in terms of the k basic moments.

First, lets define the basic moment as –

We rearrange the binomial expansion –

So we have the kth central moment given as a weighted sum of the kth simple moments –

which shows that all we need to accumulate as we walk across the sequence is the kth simple powers .

Notice the variance is now handled as a special case where k=2. Likewise k=0 corresponds to n, the element count and k=1 is the sum of elements.

c++ impl

Heres a basic impl of the above expression –

Again with usage like this –

So this way we don’t have to evaluate anything complicated in the inner loop, so its fairly quick.

You’ve probably noticed I don’t actually need the class variable sn, the sum is just the 1st order moment, and I’ve now re-factored that out too. Same for n, the element count – its just the 0th simple moment. nice.

This is good code; thanks for reminding me of it. I first saw it published in chapter 9 of Didier Besset’s “Object-oriented implementation of numerical methods”, along with a more robust but slower alternative for computing the central moments.

Thanks for the reference, haven’t read that one. My code might be similar, but was actually written from scratch independently. Im sure the approach is mentioned in other books, its probably an exercise in Knuth. In my case it was a job interview question – ‘can kth moments be calculated in one pass’ – which prompted this article.