Central moments

[LukeMathWalker]

Computation of central moments of arbitrary order.

Would it be worth to add another method that takes the mean as a parameter to avoid re-computing it if the user already has its value @jturner314?

[LukeMathWalker]

Once pairwise summation gets merged in ndarray, this will actually be a corrected pairwise two-pass algorithm for central moments.

[LukeMathWalker]

Added another method, central_moments, to compute multiple central moments up to a certain order without wasting time repeating several times the same steps (moment computation).

Useful to get an efficient kurtosis and skewness computation, given that their formulas involve central moments of different order.

[LukeMathWalker]

Added kurtosis and skewness computation.

[LukeMathWalker]

Fixed an issue in the kurtosis' test (I was comparing to SciPy's kurtosis with fisher=True instead of fisher=False).

[jturner314]

Thanks for working on this! It's especially nice that you found a numerically stable formula instead of using the naive calculation. I've created LukeMathWalker#3 with some minor changes and added a couple of comments. Everything else looks good.

[jturner314]

This epsilon is fairly large. Is it due to ndarray not having pairwise summation yet?

[LukeMathWalker]

I have tuned all epsilons to be as close as possible to 1e-12, stopping when the assertion fails.
It seems our mean estimate agrees with NumPy's with a delta of 1e-9 - I'd say that the lack of pairwise summation is a big suspect. It would be nice to have something like decimal in Rust to get a sense of the real error.

[jturner314]

I have tuned all epsilons to be as close as possible to 1e-12, stopping when the assertion fails.

Okay, thanks. That's helpful to get a better understanding of the precision.

It would be nice to have something like decimal in Rust to get a sense of the real error.

For this particular case, decimal gives the following for the standard and harmonic mean:

>>> from decimal import *
>>> getcontext().prec = 100
>>> data = [
...     Decimal('0.99889651'), Decimal('0.0150731'), Decimal('0.28492482'), Decimal('0.83819218'), Decimal('0.48413156'), Decimal('0.80710412'), Decimal('0.41762936'),
...     Decimal('0.22879429'), Decimal('0.43997224'), Decimal('0.23831807'), Decimal('0.02416466'), Decimal('0.6269962'), Decimal('0.47420614'), Decimal('0.56275487'),
...     Decimal('0.78995021'), Decimal('0.16060581'), Decimal('0.64635041'), Decimal('0.34876609'), Decimal('0.78543249'), Decimal('0.19938356'), Decimal('0.34429457'),
...     Decimal('0.88072369'), Decimal('0.17638164'), Decimal('0.60819363'), Decimal('0.250392'), Decimal('0.69912532'), Decimal('0.78855523'), Decimal('0.79140914'),
...     Decimal('0.85084218'), Decimal('0.31839879'), Decimal('0.63381769'), Decimal('0.22421048'), Decimal('0.70760302'), Decimal('0.99216018'), Decimal('0.80199153'),
...     Decimal('0.19239188'), Decimal('0.61356023'), Decimal('0.31505352'), Decimal('0.06120481'), Decimal('0.66417377'), Decimal('0.63608897'), Decimal('0.84959691'),
...     Decimal('0.43599069'), Decimal('0.77867775'), Decimal('0.88267754'), Decimal('0.83003623'), Decimal('0.67016118'), Decimal('0.67547638'), Decimal('0.65220036'),
...     Decimal('0.68043427')
... ]
>>> sum(data) / len(data)
Decimal('0.5475494054')
>>> 1 / (sum(1/x for x in data) / len(data))
Decimal('0.2179009364951495638060724591285330338229747648318613758290505954023583698639515519076072129091188061')

In Rust, we could add a dev-dependency on rug to compute the correct result to the necessary precision, but it's not necessary for this PR.

[LukeMathWalker]

That's reassuring: using the results from decimal it turns out our computation is correct up to f64::EPSILON in this specific case.

[LukeMathWalker]

We should definitely look into using something like rug to keep in check our numerical error. On the side, comparing with SciPy/NumPy reassures me that the algorithm implementation is roughly correct, even if due to the respective numerical error we do not agree on the smaller digits.

[LukeMathWalker]

I have changed the return type to match what we decided in #25 - it should be ready to be merged now @jturner314

[jturner314]

Looks good to me. The only remaining potential change is making order have type u32 or even u16 (since I can't imagine anyone wanting to calculate the 65537th order central moment). usize is fine, though.

[LukeMathWalker]

I agree - I changed it to u16 @jturner314.

[jturner314]

Now that order is u16, we can change this to let order = order as i32.

[LukeMathWalker]

True, it can't fail anymore. I'll merge once the CI has finished 👍