"This better way of computing variance goes back to a 1962 paper by B. P. Welford and is presented in Donald Knuth's Art of Computer Programming, Vol 2, page 232, 3rd edition. Although this solution has been known for decades, not enough people know about it."
I'd generalize that to "although Knuth's book has been around for decades, not enough people have read it yet."
The usual way is to define variance as the expectation of the number times its complex conjugate -- it's complex magnitude. In particular, the variance is real.
In the pseudo code, the key lines are the two after the "else" in "Push". The recurrence for the mean (first line of the two) would be unchanged.
The second line, the update for m_newS, would need a conjugate on either the first or the second parenthesized expression (in the product). Since eventually the returned variance is proportional to Sk, Sk must be real. Being real, Sk must equal its transpose. So, it must not matter which parenthesized expression you conjugate. And, you should be able to just store Sk as a real number.
To verify that it does not matter which side you put the conjugate on, just expand the defining expression for Sk (given in the article). By plugging the expression for M(k) into that for S(k), you will see the additive correction to S(k) comes out in terms of:
delta = x(k) - M(k-1)
In fact, it comes out in terms of delta times its conjugate, no matter which place you put the conjugate. (This also checks that Sk is indeed real when computed this way.) Doing this algebra is good because it allows you to simplify the complex arithmetic considerably.
I write a ton of statistical reporting code at work. We don't deal with the kind of rare cases that the author points out can cause high degrees of imprecision, but it's useful to know. The speedup alone makes changing my implementation worth it.
8 comments
[ 0.29 ms ] story [ 24.3 ms ] threadI'd generalize that to "although Knuth's book has been around for decades, not enough people have read it yet."
In the pseudo code, the key lines are the two after the "else" in "Push". The recurrence for the mean (first line of the two) would be unchanged.
The second line, the update for m_newS, would need a conjugate on either the first or the second parenthesized expression (in the product). Since eventually the returned variance is proportional to Sk, Sk must be real. Being real, Sk must equal its transpose. So, it must not matter which parenthesized expression you conjugate. And, you should be able to just store Sk as a real number.
To verify that it does not matter which side you put the conjugate on, just expand the defining expression for Sk (given in the article). By plugging the expression for M(k) into that for S(k), you will see the additive correction to S(k) comes out in terms of:
delta = x(k) - M(k-1)
In fact, it comes out in terms of delta times its conjugate, no matter which place you put the conjugate. (This also checks that Sk is indeed real when computed this way.) Doing this algebra is good because it allows you to simplify the complex arithmetic considerably.
I write a ton of statistical reporting code at work. We don't deal with the kind of rare cases that the author points out can cause high degrees of imprecision, but it's useful to know. The speedup alone makes changing my implementation worth it.
http://en.wikipedia.org/wiki/Algorithms_for_calculating_vari...