search

rust-ndarray / ndarray-stats

Statistical routines for ndarray
https://docs.rs/ndarray-stats
Apache License 2.0
201 stars 25 forks source link

Weighted var #55

Closed nilgoyette closed 5 years ago

nilgoyette commented 5 years ago

Here's the var/std version.

LukeMathWalker commented 5 years ago

Thanks for working on this! I am bit busy at the moment, I'll try get back to this during the weekend (hopefully with some answers to the points you have raised).

LukeMathWalker commented 5 years ago
  • I moved all "summary statistic" tests. Sadly github shows them as a deleted and a new file. Sorry, I should have used git mv. Ctrl-f weighted_var and weighted_std to know what I actually changed.

No worries, not a big deal.

  • I'm not too sure about ddof. I know that 0.0 and 1.0 work, but I don't actually know what others ddof do. ndarray's definition of ddof is different, so I'm probably wrong :)

I have yet to find a use for ddof outside of 0 and 1. In ndarray we have been slightly more permissive (letting it range from >0 to <n), but unless we have a concrete example I am actually ok to constrain it to [0,1].

  • How can A::from_usize(0) fails? An unwrap would be safe, no?

Well, that's exactly what we are doing :grin: We are just using expect instead of unwrap to provide a better error message in the unlikely (but still possible) event of a panic. We should use expect also when doing A::from_usize(1) and we should add this panic condition to documentation of the various new methods.

  • I picked the first 1-pass+weights algorithm that I found. I don't know if there are better or faster algorithm around, or if it accepts other ddof than 0.0 and 1.0. All I know is that it does give the same results as the 2-pass algorithm.

To notice a difference you need to evaluate it on arrays with a lot of elements, even better if they lead to somewhat badly conditioned operations. The first issue there is the summing algorithm we are using, the naive one. Using pairwise summation would greatly improve numerical stability, but there are pros and cons at the moment, considering that we do not have specialization.

Overall it looks good to me :+1: I have left a couple of comments. Could you add some benchmarks for these new methods? I am curious to see if we can make them faster using azip instead of zipping together using .iter. There are also some opportunities to use mul_add that we should probably leverage to make it faster and more precise.

nilgoyette commented 5 years ago

I added a benchmark (azip! was slower) but I'm not sure about using mul_add because criterion tells me that it's slower. It's more precise though (from 1e-10 to 1e-11). I'll add it if you want.

LukeMathWalker commented 5 years ago

That's interesting! Could you post criterion's plot? I'll try to run them later on my machine as well.

nilgoyette commented 5 years ago

I replace

mean += (w / weight_sum) * x_minus_mean;
s += w * x_minus_mean * (x - mean);

with

mean = (w / weight_sum).mul_add(x_minus_mean, mean);
s = (w * x_minus_mean).mul_add(x - mean, s);

cargo bench gives something like

   /10 Change within noise threshold.
  /100 +2.6839%
 /1000 +3.0680%
/10000 +1.9234%
LukeMathWalker commented 5 years ago

The performance impact looks even more severe on my machine: image

I'd say we can go without mul_add for now and then revise if we find out we need to squeeze more precision.

LukeMathWalker commented 5 years ago

Same story for Zip (quite interesting actually): image

I'd say we are good to go :muscle: Thanks for working on this!