ENH: add `quantile`

[lucyleeow]

It was decided not to add quantile to array-api (data-apis/array-api#795), partly because scipy added array API support to their quantile (scipy/scipy#22352) and partly due to complications (many methods, weight support, nan support).

This RFC is mostly motivated by usage in scikit-learn (as seen from data-apis/array-api#795 (comment) its not used much in other packages).

Although scipy does have an array API quantile, due to our min scipy version, we are currently not able to use it (e.g., scikit-learn/scikit-learn#31671).

More generally, it seems to be a commonly enough used function (which is implemented many array libraries) to warrant considering adding it to here, and having a consistent, single source.

Not sure who would be interested but maybe:
cc @betatim @ogrisel @lorentzenchr (and @rgommers who suggested it)

[lucascolley]

PR welcome! We can always deprecate or declare legacy once the SciPy implementation is sufficiently old enough

[betatim]

What do you think of basically copying the scipy implementation, checking what version of scipy is installed and then based on that decide to use scipy or the version in array-api-extra?

[lucascolley]

yes, the SciPy version should work for https://github.com/data-apis/array-api-extra/blob/main/src/array_api_extra/_lib/_funcs.py (pending any test failures), and delegation to existing backends (including SciPy in this case) can happen in https://github.com/data-apis/array-api-extra/blob/main/src/array_api_extra/_delegation.py.

Contributor docs at https://data-apis.org/array-api-extra/contributing.html !

[lorentzenchr]

The version in scipy,stats.quantile does not (yet?) support (sample) weights, while scikit-learn uses them in several places (trend: growing, see _weighted_percentile). The real benefit would be to have an Array-API compatible version with weights. (On the other hand, 2 methods would do for me: inverted_cdf, linear).

[lucyleeow]

Thanks all, I didn't expect such a quick positive response! And thanks for tackling this @betatim !

@OmarManzoor raised a good point about most array libraries already implementing their own quantile (which may be faster than our implementation). However, there are many ways to calculate quantile.

The following packages implement all 9 methods from the Hyndman and Fan paper ¹ :

• numpy.quantile
• xarray.DataArray.quantile

These packages implement the older, mostly discontinuous ('linear', 'lower', 'higher', 'midpoint', 'nearest') methods that numpy used to have (but still support for backwards compatibility), which are not described in the Hyndman and Fan paper (except 'linear'):

• torch.quantile
• cupy.quantile
• jax.numpy.quantile
• dask.array.quantile
• tfp.stats.quantiles

In light of this, and for future proofing, having our own implementation would be great!

I am also hopeful that implementing here means that we have a single path towards pushing for weight and nan support 😬 both of which we would like to use in scikit-learn.

Numpy added weight support for "inverted_cdf" in numpy/numpy#24254 (thanks @lorentzenchr ). Ideally we would add support for another method, whose output matches that of median (among other desirable properties). "linear" (continuous) and "average_inverted_cdf" (discrete) have both been suggested (#340 (comment), numpy/numpy#24254 (comment))

NaN support is possibly less compelling (but it would be nice). We've recently added nan support to SplineTransformer but would like to do so for KBinsDiscretizer and maybe RobustScaler and QuantileTransformer.

Footnotes

1. R. J. Hyndman and Y. Fan, “Sample quantiles in statistical packages,” The American Statistician, 50(4), pp. 361-365, 1996 ↩

[betatim]

I've not educated myself about the (many) different ways of calculating quantiles. I think two things are desirable: (1) use "native" implementations as much as possible and (2) provide a quantile function that supports what users want. It sounds like there is some tension between the two :-/

Things we should IMHO not aim for: (1) provide as many different methods for computing quantiles (if it happens for free then that is fine, but I'm not sure it is worth putting effort into this) and (2) changing what quantile does once it exists in array-api-extra (let's try and "get it right" from the start so that we don't have to deal with deprecations and all that)

In particular sample weights seem to not have a lot of support in libraries. Does this mean we should not rely on the underlying array libs and "roll our own"?

[lucascolley]

There has been discussion in SciPy about how there may be various competing interpretations for a weights parameter, hence why stats.quantile does not have one (yet). I can try to find the discussions if that is helpful.

[lucyleeow]

In particular sample weights seem to not have a lot of support in libraries.

Sample weights is no simple feat. Just see this PR in numpy, trying to add sample weights, which was active for almost 6 years (!!) and then closed, unmerged.
Sebastian outlines the problems well, here.

When you see the terms fweights, aweights etc, I think they are referring to the stata terms (who coined 'analytical' weights). It does not help that people call them different terms, here is a summary in case it is useful:

• frequency weights
• analytical (stata) weights, aka precision weights, aka reliability weights, aka inverse variance weights (ref)
• probability weights, aka sampling weights

I have (suprisingly) found a decent source describing these different interpretations of weight and how they are handled: http://www.stata.com/manuals13/u20.pdf#u20.23Weightedestimation (see section 20.23) It talks about weights generally, and not specifically wrt quantile.

My limited understanding is that the only difference in handling is whether the weights are normalized to n or not. This is mentioned in R Hmisc (I feel like we can trust R to get stats right, right?): https://search.r-project.org/CRAN/refmans/Hmisc/html/wtd.stats.html (see normwt parameter). Though honestly, I'd need to look more deeply into this.

use "native" implementations as much as possible

I had the same thought initially. We could even get past the different implementations part because I think they all implement 'linear', but because we sort of 'need' weight and nan support, I think we have to roll our own?

[betatim]

What about using the inverted_cdf method that is implemented in numpy and supports weights? That way we have a good implementation to base ours on, we get sample weight support. Is the inverted_cdf method "too weird"/why is it not the default in Numpy (except for backwards compatibility issues)?