ENH: improve siegelslopes via pythran

[charlotte12l]

The Pythran version is 10x faster than the original Python version. cc @rgommers @serge-sans-paille

Timing codes:

import numpy as np
from timeit import timeit
import scipy.stats.mstats as mstats

x = 2 * np.arange(100)
y = 5 * x - 3.0

f1 = lambda : mstats.siegelslopes(x,y)
print(timeit(f1, number=100)) #0.020322154000000037 vs 0.271933536 10x faster

[rgommers]

Thanks @charlotte12l! This looks straightforward, except for the same discussion as in gh-14429 about the many signatures. Let's resolve that on gh-14429 and then see if we're happy with this PR as it is now.

[charlotte12l]

Thanks @charlotte12l! This looks straightforward, except for the same discussion as in gh-14429 about the many signatures. Let's resolve that on gh-14429 and then see if we're happy with this PR as it is now.

I just left the public function and doc in Python, as we discussed before.

[rgommers]

Thanks @charlotte12l, this is looking quite good now. One more design issue to deal with about dtypes.

[rgommers]

The one thing to look at still is this type annotation. The types are not constrained in the public Python function, so this won't work for, e.g., float32 or int16. I think we still need to figure out exactly how to deal with dtypes in these cases.

[rgommers]

There's a balance to strike between binary size, performance, and maintainability.

[mdhaber]

It doesn't look like the original Python function had support for different data types in mind, as x is explicitly cast to float at the beginning. Were we hoping to fix this in the Pythran version?

(Just to confirm, since you don't mention complex dtypes here - contributors probably don't often have complex generalizations in mind when they write stats functions. Do we have a policy about support for complex dtypes in stats functions?)

For now, I think that a Pythranized version of this function would be useful even if it were just for float64 (and maybe also float32). The benefit of integer versions seems less compelling here. There could be some memory savings for int8 and int16, but it will probably be slower as the type needs to be promoted to do the division.

If we want to be more careful about dtypes throughout stats in the future, we could come back to it. In the interest of seeing these PRs through, would supporting just those one or two types make sense?

[mdhaber]

This looked close, so I went ahead and exported only float32 and float64, and I use np.result_dtype to determine which of those dtypes to convert to before calling the function. Does that look ok @rgommers?

[mdhaber]

This PR renames _hypotests_pythran to _stats_pythran (since it's not just for hypothesis tests) and copy-pastes (with mild modifications) the innards of siegelslopes into _stats_pythran for a nice speed boost It seemed to be held up by a discussion of what data types it should be compiled with; I went ahead with float32 and float64 (#14430 (comment)). I'll go ahead and merge, but please let me know if I missed anything @rgommers @serge-sans-paille and I can fix it in a new PR.