extra::stats suffers from catastrophic cancellation

[huonw]

It is using the naive algorithm which is numerically unstable; it should be using something like Kahan summation.

extern mod extra;
use extra::stats::Stats;

fn main() {
    // prints 0, but should be 1
    println!("{}", [1e20, 1.0, -1e20].sum());
}

[huonw]

I'm happy to mentor someone through this as a relatively easy first bug. (Feel free to talk to me here, or ping me (nick: dbaupp) on IRC.)

[g3xzh]

If I've understood correctly, I have to implement Kahan summation here.
Probably to change this part |p,q| p+q*.
Right?

[huonw]

Yes, that's exactly correct; you can either convert the .fold into a for loop or just make the fold pass through a tuple with the two piece of state, i.e. .fold((0.0, 0.0), |(sum, c), q| { ... (new_sum, new_c) }).

Also, another neat thing to do would be to make the variance calculation more numerically stable, e.g. using this or this.

[huonw]

(Also, sorry for repeatedly missing you on IRC; I'm in Australia which means that our times seem to not line up. :( )

[g3xzh]

@huonw
Regarding the IRC, that's fine. :)
Anyway, if you don't mind, I take a shot and work on it this week.

[huonw]

I certainly don't mind, go for it! Feel free to attempt to catch me in IRC if you have any questions (or just throw them here, or email me (should be visible in my github profile)).

[g3xzh]

Hey @huonw,

I have implemented the compensated summation but it has not done any impact.
When [1e20, 1.0, -1e20].sum() the 1.0 is too small relative to 1e20, making it negligible, thus 1e20+1.0=1e20.
Another way to fix it, is to sort the vector by its elements' absolute values and then do the summation.
However, the sort will add significant overhead and I am not sure if it is efficient and worthwhile.

[huonw]

Hm, Python has math.fsum which works perfectly:

>>> import math
>>> math.fsum([1e20, 1, -1e20])
1.0

so we may need to use a different algorithm.

[g3xzh]

OK, I will take a look at this :D.
BTW: that's my implementation.
http://pastebin.com/QqjVWLR1

[huonw]

The implementation of fsum, which appears to be the msum algorithm from my previous link.

It would be good for the documentation on sum to mention (a) that this method sacrifices performance at the altar of accuracy, and (b) to cite the paper that it apparently comes from (although I have not read it).

(That implementation of Kahan summation looks good.)

[huonw]

Fixed by #10927.