issue_comments

Merging now, well done everyone!

Turns out align does not entirely do what I thought it does. For e.g. da1 * da2 there is an align call for the DataArray and then broadcasting happens on the level of the Variable - which we skip here. Therefore, I updated the PR (5c251e0) - to use xr.broadcast to ensure da and weights have the same shape.

Notes 1. The tests should fail because I updated the skipna logic and this was also "wrong" (inconsistent with the rest of xarray) in xr.Dataset.quantile - thus this needs #6303 first. 2. broadcast does an outer join which is kind of pointless here because it adds NaN values that are then removed again (#6304)

Thanks, for checking - I didn't anticipate the missing core dim... Seems the align may be required, but it should be possible to leave out the shape check. I'll have a closer look (also how dot does it - because this is what's used for the other weighted methods).

The alignment of data and weights is not done automatically. So I agree this would be ideal, but I'll need some guidance to make it happen without the align call.

This is the error I get python ValueError: operand to apply_ufunc has required core dimensions ['dim_0', 'dim_1'], but some of these dimensions are absent on an input variable: ['dim_1']

@mathause I'm not 100% confident the methods other than linear work as expected, so I suggest we do not expose method until a robust test suite can confirm results are valid.

I have added a test to verify that using equal weights with the different interpolation methods that this PR supports would work (i.e. would yield the same results as np.quantile, with the corresponding methods). It is skipped however because the method argument is not currently exposed in the API (it would be in future work, ideally).

We have local changes that expose method in the public API, and compare results against numpy for the case with uniform weights to confirm results are identical.

Do you want us to push those changes here ?

If you are confident they work with the weights... numpy now recommends using median_unbiased or normal_unbiased, so would be nice if we have them.

We have local changes that expose method in the public API, and compare results against numpy for the case with uniform weights to confirm results are identical.

Do you want us to push those changes here ?

@mathause This PR goes beyond what is currently implemented in numpy. For now, all weighted quantiles PR on numpy are more or less based on "linear" method (method 7) and none have been merged. I plan to work on integrating weights with the other interpolation methods but don't have the time right now. I'll probably pick some ideas from here.

As for the numerics here, everything looks good to me. The only limitations I can see are: - This only handles sampling weights, which is fine I guess. - Some interpolation methods are missing, they can be added later. - ~A nan_weighted_quantile could also be interesting to add~

Correct, it's not exposed yet because I don't have the bandwidth to create tests for all the different methods.

The equal weight case could be tested against numpy fairly easily though.

Good point, yes the keyword should be named method as well for consistency (but as far as I can see it is currently not exposed?).

Sorry for the slow pace - this PR merits another review. I can probably look at the code but I am not sure about the numerics. I think it's important to align with numpy but their code is not yet merged (AFAIK). Maybe @bzah knows if this would be consistent with numpys approach?

@mathause Should this PR be amended to account for #6108 ?

FYI, weighted quantiles topic will be discussed in numpy's triage meeting of today (17:00 UTC). I'm not a maintainer but I'm sure you are welcomed to join in if you are interested. Meeting information: https://hackmd.io/68i_JvOYQfy9ERiHgXMPvg

We discussed this in the Triage meeting a bit, and settled on "we want this". We are OK with the sorting approach and limiting e.g. only to aweights, that may need settling 100%, but generally using aweights (and maybe fweights) sems good.

EDIT: Oops, meant to post on the NumPy PR, but OK here as well, I guess ;)

I had the same interrogation. My guess is that the DataArray.quantile and DataArrayWeighted.quantile should have similar signatures.

The way I see it is that another PR would bring the DataArray.quantile signature up to date with numpy 1.22, which includes the 9 h types, and harmonize DataArrayWeighted.quantile with it.

@huard's latest commit modifies the algorithm so that it uses Kish's effective sample size, as described in the blog where the algorithm comes from: https://aakinshin.net/posts/kish-ess-weighted-quantiles/, which seems to solve the problem mentioned by @mathause.

Also he adds support for the interpolation types 4 to 9 (those that share a common way of computing Qp, as described here: https://en.wikipedia.org/wiki/Quantile#Estimating_quantiles_from_a_sample

Many thanks for calling this to my attention (via #5870)! Following.

@mathause About this:

I did some tries and got an unexpected result:

```python data = xr.DataArray([0, 1, 2, 3]) weights = xr.DataArray([1, 0, 1, 0]) data.weighted(weights).quantile([0.75])

np.quantile([0, 2], 0.75) ```

Can you double-check? Or do I misunderstand something?

My latest commit should fix (and test) this.

@mathause Thanks for the many excellent suggestions! After having removed the for loop the way you suggested, I tried to address this:

The algorithm is quite clever but it multiplies all elements (except 2) with 0 - this could maybe be sped up by only using the relevant elements.

At first I thought that something like this could work:

python w = np.diff(v) nz = np.nonzero(w) d2 = np.tile(data, (n, 1)) r = w[nz] * d2[nz] r = r[::2] + r[1::2]

The problem however is that it turns out that w's rows sometimes have one element only, instead of two (which is when an h coincides exactly with a weight value, instead of lying between two). Given that difficulty, my impression is that it's not really solvable, or at least not in a way that would result in a more efficient version.