ENH: Port null-distribution special functions from scipy.stats to xsf

[fbourgey]

Following @mdhaber's comment and @steppi's suggestion.

Motivation

Several null distributions of test statistics are currently implemented in scipy.stats, even though they are essentially special functions. The issue tracks and coordinates porting these functions to xsf.

Moving these functions to xsf would:

• make them more reusable by other libraries
• leverage Special’s infrastructure (delegation, lazy_apply, etc.)
• enable efficient compiled (C/C++) implementations
• improve performance for repeated CDF/SF evaluations at many statistic values for the same sample size, where Python-level caching is often insufficient

Proposed starting point

• _cdf_cvm_inf: asymptotic (infinite-sample-size) null distribution CDF of the one-sample Cramér–von Mises statistic
  • ENH: Implement cvm_cdf_inf from scipy.stats._hypotests #97

This would serve as a pilot for additional ports.

Other asymptotic or exact null distributions that could follow

• Kolmogorov-Smirnov one-finite-sample, two-sided distribution
• Cramér–von Mises one-sample
  • ENH: Implement cdf_cvm from scipy.stats._hypotests #104
• Cramér–von Mises two-sample
  • ENH: Implement the exact p-value of the Cramér-von Mises two-sample test #108
  • ENH: Add exact Cramér-von Mises two-sample p-value gufunc #118
• Ansari
• Wilcoxon
  • ENH: Add CDF and SF of the Wilcoxon distribution #124
• Mann–Whitney U
• Kendall’s tau
  • ENH: Implement Kendall discordance #110
• Kolmogorov–Smirnov two-sample exact
• von MIses distribution
  • ENH: Implement CDF of the von Mises distribution #111
• Poisson binomial distribution
  • ENH: implement poisson_binom from scipy.stats._stats_pythran #101
  • ENH: Add poisson_binom_all gufunc kernel #113
  • ENH: Add addional gufunc kernels for poisson_binom distribution #117
• Bivariate normal distribution
  • ENH: Implement CDF of the bivariate normal distribution #112
• Trivariate normal distribution
• Multivariate normal distribution (also _quato)
• Bivariate t distribution
• Trivariate t distribution
• Multivariate t distribution
• Noncentral hypergeometric distributions

Sample statistics:

These are really tough (maybe impossible) to express efficiently in terms of the array API.

• Kendall's Tau

[fbourgey]

@mdhaber @steppi feel free to update/correct the different functions.

[mdhaber]

Thanks for opening this @fbourgey! I've updated the list with all I could find right now. A few questions:

1. Can you or @steppi outline how these can behave like ufuncs but cache results (at least within a given call to the ufunc)?
2. Some of the ones I've added above, like the multivariate t and normal distributions, require randomness. It is important to be able to use a NumPy Generator (or its lower-level interface) to control it. How can that be handled?
3. There are some sample statistics - the statistics themselves, not the null distributions - that are difficult to express in terms of the array API without tight loops. Is there room for that sort of thing in xsf?

[fbourgey]

Thanks, @mdhaber, for updating the list. I'll defer to @steppi on those questions. Once we have a clearer picture of how the conversion to xsf should be handled, based on the example #97, I'll start working on the others.

[steppi]

1. Can you or @steppi outline how these can behave like ufuncs but cache results (at least within a given call to the ufunc)?

It can be done with custom ufunc loops. Usually loops just loop over the input arrays and compute values elementwise and populate the output array, but you can run pretty much arbitrary code C/C+= code in them. We could implement a cached_ufunc implementation in xsf/numpy.h with a custom loop in place of the typical style loop that's used here:

https://github.com/scipy/xsf/blob/main/include/xsf/numpy.h#L790-L813

There would be a need to think about how to write a sane interface for instantiating such cached ufuncs. I had once discussed custom ufunc loops with Irwin, and he thought some of my more complex ideas were too much trouble to bother with, but he seemed to think adding a cache would probably be fine.

Getting this right wouldn't exactly be a simple undertaking. It's certainly doable though.

2. Some of the ones I've added above, like the multivariate t and normal distributions, require randomness. It is important to be able to use a NumPy Generator (or its lower-level interface) to control it. How can that be handled?

do you need to use the generator internally in the kernel or is it fine to pass in an array of randomly sampled values? The latter is doable, the former would be a non-starter in xsf I think.

3. There are some sample statistics - the statistics themselves, not the null distributions - that are difficult to express in terms of the array API without tight loops. Is there room for that sort of thing in xsf?

I have no qualms about adding them.

[mdhaber]

1. Whatever the implementer and @steppi agree on for starters is fine with me. OK with me to start without a cache and add later, if that's easier.
2. I remember from the initial implementation of MVN/MVT, it would be possible to pass in uniform random data. There is a rejection component, but I think we had a plan for dealing with that. It's OK to leave those at the bottom of the list.
3. Great. I'll add a list.

[steppi]

1. Whatever the implementer and @steppi agree on for starters is fine with me. OK with me to start without a cache and add later, if that's easier.

I should have pointed out that the scalar kernels one would write if not using caching wouldn't have to be changed to use caching through the ufunc loops. So if @fbourgey writes them without using caching, I can go back later and add the custom loops that allow caching and set things up.

2. I remember from the initial implementation of MVN/MVT, it would be possible to pass in uniform random data. There is a rejection component, but I think we had a plan for dealing with that. It's OK to leave those at the bottom of the list.

Cool. It seems like the rejection component could be handled by passing in more random points than needed by a comfortable enough margin.

[fbourgey]

For KS, do we port all methods (pdf/cdf/sf/ppf/isf)?

[mdhaber]

I guess I only need the ones that are used in the test. Up to Albert whether he wants the rest, but maybe we at least start with the cdfs/sfs for all the distributions, and come back to other methods later if needed.

[fbourgey]

I guess I only need the ones that are used in the test. Up to Albert whether he wants the rest, but maybe we at least start with the cdfs/sfs for all the distributions, and come back to other methods later if needed.

It seems there is already a kolmogorov.h file. This was added by @steppi in 2024 but I think the implementation differs from SciPy. @steppi can you confirm and should this be reimplemented?

Should we follow the same pattern and translate the whole class to C++? I was thinking of starting with Ansari in a new file ansari.h

[mdhaber]

My suggestion was to start with the methods actually used by the hypothesis tests. If they are all done for KS tests, then it would be great to move onto Ansari. Coming back later for completeness is an option, but my initial interest is just getting the null distributions used in the hypothesis tests into XSF.