Repository navigation
ENH: add quantile #340
Description
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)
Activity
lucascolley commentedon Jun 30, 2025
PR welcome! We can always deprecate or declare legacy once the SciPy implementation is sufficiently old enough
betatim commentedon Jun 30, 2025
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 commentedon Jun 30, 2025
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 !
[-]RFC Add `quantile` ?[/-][+]ENH: add `quantile`[/+]lorentzenchr commentedon Jun 30, 2025
The version in
scipy,stats.quantiledoes 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 commentedon Jul 2, 2025
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 1 :
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'):
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
SplineTransformerbut would like to do so forKBinsDiscretizerand maybeRobustScalerandQuantileTransformer.Footnotes
R. J. Hyndman and Y. Fan, “Sample quantiles in statistical packages,” The American Statistician, 50(4), pp. 361-365, 1996 ↩
betatim commentedon Jul 2, 2025
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
quantiledoes 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 commentedon Jul 2, 2025
There has been discussion in SciPy about how there may be various competing interpretations for a
weightsparameter, hence whystats.quantiledoes not have one (yet). I can try to find the discussions if that is helpful.lucyleeow commentedon Jul 2, 2025
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,aweightsetc, 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: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
normwtparameter). Though honestly, I'd need to look more deeply into this.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 commentedon Jul 2, 2025
What about using the
inverted_cdfmethod 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 theinverted_cdfmethod "too weird"/why is it not the default in Numpy (except for backwards compatibility issues)?lucyleeow commentedon Jul 3, 2025
Honestly, for large
nI feel like any method would be okay (and may be considered within statistical uncertainty). See scikit-learn/scikit-learn#17370 (comment) and scikit-learn/scikit-learn#17370 (comment). With that said:With the general equation$Q_{i}(p)=(1-\gamma)X_{(j)}+\gamma X_{(j+1)}$
inverted_cdfis probably the simplest method:Index is defined as$j\le pn\lt j+1$ ($n$ is sample size and $p$ is percentile)
then:
where$g=pn-j$ (essentially the 'fractional' component above $pn$ )
But it does not have any nice symmetry properties (which is why we added
_average_weighted_percentilein sklearn, see scikit-learn/scikit-learn#17370 (comment)). In particular:(Honestly I don't intuitively have a feel on how these are different, except I guess for continuous quantile methods, the cdf/plotting position formula would vary, so they don't mean the same?)
Also, Q(0.5) is only equal to sample median when n is odd.
average_inverted_cdfsatisfies the first symmetry property and Q(0.5) is always equal to sample median.linearsatisfies the second symmetry property but not the first and Q(0.5) is always equal to sample median. It's also continuous.Practically - weight support for
inverted_cdfwould probably be good enough (🤷) but weight support foraverage_inverted_cdf, orlinearwould probably be nicer and help some of our tests pass (e.g., scikit-learn/scikit-learn#29818 and scikit-learn/scikit-learn#16298) 😅edit: actually, when
sample_weight=Nonewe usenp.percentilewhich defaults to 'linear' method, so it would be nice if the calculation with weights didn't differ too much from the calculation without weights.lucyleeow commentedon Jul 10, 2025
@lucascolley did you manage to find discussions in SciPy around weight support? Maybe there are people more knowledgeable in stats in SciPy? (ref #340 (comment))
12 remaining items
lucyleeow commentedon Oct 21, 2025
Not sure if delegating is possible, it relies on all packages having consistent behaviour, which we'd need to check. Just taking a quick look at tensorflow tfp.stats.quantiles - and they take
num_quantiles, instead ofq.Also it seems that scipy and numpy
axisbehaviour also differs: scipy/scipy#23832cakedev0 commentedon Oct 21, 2025
I think tensorflow is not a supported backend of
array-api-extra(see https://github.com/data-apis/array-api-extra/blob/main/src/array_api_extra/_lib/_backends.py).I already wrote the delegation code for numpy, cupy, torch, dask and jax, and based on the docs, it should be fine (but we'll see if the tests pass).
Ah, thanks for the pointer, I've been scratching my head on this broadcasting thing for a while, this should help ^^
lucyleeow commentedon Oct 21, 2025
Ah interesting, I'm sure they had a good reason. Tensorflow is listed as a stakeholder (in array API standard): https://data-apis.org/array-api/latest/purpose_and_scope.html#stakeholders
If we aren't interested in tensorflow, then seems reasonable, depending on how much tweaking is required for each individual package (as this may be difficult to maintain going forward). Do you want to open a draft PR and we can take a look?
lucyleeow commentedon Oct 21, 2025
Out of interest, why exclude tensorflow @lucascolley ?
lucascolley commentedon Oct 21, 2025
I don't think we are concerned with supporting tensorflow in array-api-extra at least, since it isn't supported by array-api-compat: https://data-apis.org/array-api-compat/supported-array-libraries.html. I'm missing historical context for the relationship between tensorflow and the consortium.
lucyleeow commentedon Oct 21, 2025
For a general guide to broadcasting, this is nice: https://quansight-labs.github.io/ndindex/indexing-guide/other-topics.html#broadcasting
cakedev0 commentedon Oct 21, 2025
Sorry, what I meant by "broadcasting thing" was not broadcasting in general but what was done in numpy/scipy for the quantile function related to broadcasting.
In numpy, there is an error:
ValueError: q must be a scalar or 1d. (but it's a bit weird, because it works with 2D q, see issue above).So I propose to support only scalar and 1d for
q. This should be enough for what is needed in scikit-learn right? And it would make the implementation quite easier I think.lucascolley commentedon Oct 21, 2025
yes, a good rule of thumb is to only implement a sensible subset of the APIs based on what we need, for more complex functions.
cakedev0 commentedon Oct 21, 2025
Ok great! I should be able to make my draft PR readable and tested soon 👍
cakedev0 commentedon Oct 22, 2025
I'm extending the scope of my PR. In the end, I prefer doing most of it in one-go to reduce context switchs (and I feel this will be true for reviewers too), because it's not so hard now that I have a clear understanding of the scope and the desired shapes/broadcasting behavior.
New scope:
"linear", "inverted_cdf", "averaged_inverted_cdf"nan_policy="omit"for weighted quantile. And usual behavior isnan_policy="propagte"(numpy's behavior for instance). So I might implement both. Or just always omit. What do you think?lucyleeow commentedon Oct 23, 2025
A thought (which may or may not be much of a concern): will it be odd that we only support 'linear' if the array library is one of the supported backends but 'average_inverted_cdf' (and not 'linear') otherwise?
I does mean when writing code, I probably wouldn't use 'linear' because it's not as widely supported as 'average_inverted_cdf'.
We could of course add support for 'linear', but in that case, we may as well add support for the other continuous methods as well.
But I am not 100% we have full agreement that limited support of method(s) would be the way to go, WDYT @betatim @ogrisel and maybe @lesteve
cakedev0 commentedon Oct 23, 2025
I'm not sure to understand what you mean. For me, the case "otherwise" is empty, because in my understanding "otherwise" means no
array-api-compatand hence no ways to get a workingxp.lucyleeow commentedon Oct 23, 2025
Ah I may have some mis-understanding. In array API extra, are the 'extra' functions ONLY for array libraries supported by array-api-compat?
I thought they were meant to be written general enough to support 'any' library that supports the standard?
cakedev0 commentedon Oct 23, 2025
Hum... My understanding is that any library that supports the standard can be trivially added to
array-api-compatso I think it's fair to assume that the set "libraries that support the standard but aren't inarray-api-compat" is empty.But that's just my understanding, so I'll let @lucascolley or someone else confirm.
cakedev0 commentedon Oct 23, 2025
They are. They are tested against
array_api_strictandarray_api_strictest(unless explicitly skipped, a few tests skip strictest, but none skip strict I think).