✔️ Linting Passed

All linting checks passed. Your pull request is in excellent shape! ☀️

_{Generated for commit: f0e999e. Link to the linter CI: here}

I think this is ready for review, maybe @ogrisel @betatim ?

I looked into the failure from #29641 (comment) (I decided it was more relevant to detail here, and better that it does not get lost in that much bigger PR).

tl;dr the difference is due to floating point precision error, and the new implementation more closely matches how quantiles are implemented in numpy.

Using the test from #29641 (comment), I made a 'minimal' (enough) reproducer. Note that the percentile indices [9, 26, 36] differ between the 2 methods, but the example below just shows 9.

def test_weighted_percentiles():
    rng = np.random.RandomState(63)

    n_samples = 200
    X = rng.randn(n_samples, 3)
    sw = rng.randint(0, 5, size=n_samples)

    # repeat and compute using numpy non-weighted
    X_repeated = np.repeat(X, sw, axis=0)
    col_data = X_repeated[:, 0]
    sort_idx = np.argsort(col_data)
    col_data = col_data[sort_idx]
    percentiles = np.linspace(0, 100, num=49 + 1)
    percentiles = percentiles[1:-1]

    np_percentile = np.percentile( percentiles[9], method="averaged_inverted_cdf")
    print(f'{np_percentile=}')

    # Do not repeat and compute using our weighted_percentiles with weights
    col_data = X[:, 0]
    nnz_sw = sw != 0
    col_data = col_data[nnz_sw]
    sw = sw[nnz_sw]
    sort_idx = np.argsort(col_data)
    col_data = col_data[sort_idx]
    sw = sw[sort_idx]

    old = _averaged_weighted_percentile(col_data, sw, percentiles[9])
    print(f'{old=}')

    new = _weighted_percentile(col_data, sw, percentiles[9], average=True)
    print(f'{new=}')

In this case ([9]), the 'new' method (this PR): adjusted percentile rank (i.e. quantile * cumulative sum) ends up 80.0 exactly (I checked and it matches to 15 decimal places). This matches a value in the weighted_cdf (note sample weights are all integers), so fraction above is 0, and we take the average of the 2 adjacent values.

With the old method, when we take _weighted_percentile of 100-percentile, adjusted percentile ends up being 312.00000000000006. Technically it should be 312 (i.e. cumulative sum - adjusted percentile rank of percentile (80)). This extra bit just pushes us to the next index, so _weighted_percentile(percentile) and _weighted_percentile(100-percentile) end up being the same number - but 100-percentile 'should' have been the adjacent value (assuming the forward adjusted rank value was slightly more accurate).

The other indices, [26, 36], have a similar problem where the adjusted percentile of percentile and 100-percentile don't add up exactly to the cumulative sum (weighted_cdf[-1]), due to floating point precision, and this also just happens to cause searchsorted to pick the adjacent index, just because of the values that were in the weighted_cdf.
I guess the extra calculations (100-percentile then the adjusted_percentile_rank = percentile_rank / 100 * weight_cdf[..., -1]) gives more chance for error.

It's tricky to compare to numpy because they don't have a weighted implementation, and without weights you can just use n to calculate - which is what they do (ref):

_QuantileMethods = {
    # --- HYNDMAN and FAN METHODS
    # Discrete methods
    'inverted_cdf': {
        'get_virtual_index': lambda n, quantiles: _inverted_cdf(n, quantiles),  # noqa: PLW0108
        'fix_gamma': None,  # should never be called
    },
    'averaged_inverted_cdf': {
        'get_virtual_index': lambda n, quantiles: (n * quantiles) - 1,
        'fix_gamma': lambda gamma, _: _get_gamma_mask(
            shape=gamma.shape,
            default_value=1.,
            conditioned_value=0.5,
            where=gamma == 0),
    },

But ultimately I think this calculation (n * quantiles) - 1 matches the new implementation more than reversing the array and calculating the 100-percentile.

I'll fix the merge conflict and address review now, thanks for review!

Thanks @ogrisel ! New commit only adds comments

Thanks for the review again @ogrisel , it reads much better now.

@ogrisel @virchan I've added the constant multiplier test and re-ran the CUDA CI.

@ogrisel do you think we should ping someone else for review?

Let me trigger the CI with all random seeds on the new test and then merge if green.

@lucyleeow the new test fails for some random seeds. It might have discovered a real bug?

Let me trigger the CI with all random seeds on the new test

Good call!

The test failures are the new constant multiplier test with float, due to cumulative sum error, which I did suspect:

it certainly should be true for integer c but for float c, it is possible that floating point precision problems may cause inequality

We've also seen the same problem here: #30787 (comment)

I double checked one case - [93-20-True-0.3] and can confirm.
For the raw percentile calculation adjusted_percentile_rank (percentile_rank / 100 * weight_cdf[..., -1]) is 7.0, falling exactly on a data point.
For the 0.3 multiplied calculation adjusted_percentile_rank is 2.0999999999999996, falling just below 2.1, making is_fraction_above to be True (instead of False, which it was in the first case).

Looking at the failures, it does pass most of the time, which is nice.

I will change the test to check 2 integer multipliers instead.

I've re-triggered the all random seeds

@ogrisel are you happy for this to go in?

@ogrisel gentle ping

Sorry for the slow feedback, thanks very much!