DOC: specify behavior of `_weighted_percentile` with `average=True`

### Describe the bug

When using `sklearn.utils.stats._weighted_percentile(..., average=True)` with floating-point weights and such that the "average code path" is expected to be taken (i.e. [those lines](https://github.com/scikit-learn/scikit-learn/blob/aae146987f3ded82dc631dc1f194863917c575b4/sklearn/utils/stats.py#L181-L185)), the result in unpredictable:
- sometimes you do get the average
- sometimes you get the value just below
- sometimes you get the value just above

I think this is not easy to fix: floats have rounding errors and `average=True` has to rely on exact equality.

Note that the functions behave properly for intergers (as long as weights.sum() < 2^52 or something like that). Indeed, even if integers are converted to floats in the function, the 52 bits are enough exactly represent intergers up to ~2^52 which is more than enough.

So maybe one way to fix this would be to transform/scale the weights in a way that summing them causes no rounding errors. This probably means that all weights have to share the same exponent part in their floating-point representation, and that summing them don't overflow the mantis part. I will try to explore this idea.

### Steps/Code to Reproduce

Imports:
```Python
from sklearn.utils.stats import _weighted_percentile
import numpy as np
```

Three small snippets:
### 1
```Python
x = np.arange(10)
w = np.ones_like(x)
print(
    _weighted_percentile(x, w, average=True),
    _weighted_percentile(x, w * 3 * np.pi, average=True),
    _weighted_percentile(x, w * 30 * np.pi, average=True)
)
```
prints 4.5 4.0 5.0 => you would expect 4.5 4.5 4.5

### 2

```Python
x = np.arange(10)

for _ in range(100):
    w_half = np.random.randint(1_000_000, size=x.size // 2)
    w = np.concat((w_half, w_half))
    # will never fail:
    assert _weighted_percentile(x, w, average=True) == 4.5

for _ in range(100):
    w_half = np.random.rand(x.size // 2)
    w = np.concat((w_half, w_half))
    # will sometimes fail:
    assert _weighted_percentile(x, w, average=True) == 4.5
```

### 3

```Python
x = np.arange(100_000)
w_half = np.random.rand(x.size // 2)
w = np.concat((w_half, w_half))
_weighted_percentile(x, w, average=True)
```
=> (almost?) never outputs 49999.5 as expected, but always 49999.0 or 50000.0
The rounding error is more likely in this example than with `x.size == 10`.

This shows that basically, **for large arrays, the option `average=True` and `average=False` are basically equivalent**.

### Expected Results

See comments for each snippet

### Actual Results

See comments for each snippet

### Versions

```shell
System:
    python: 3.12.11 (main, Aug 18 2025, 19:19:11) [Clang 20.1.4 ]
executable: /home/arthur/dev-perso/scikit-learn/sklearn-env/bin/python
   machine: Linux-6.14.0-32-generic-x86_64-with-glibc2.39

Python dependencies:
      sklearn: 1.8.dev0
          pip: None
   setuptools: 80.9.0
        numpy: 2.3.3
        scipy: 1.16.2
       Cython: 3.1.4
       pandas: 2.3.2
   matplotlib: 3.10.6
       joblib: 1.5.2
threadpoolctl: 3.6.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 16
         prefix: libscipy_openblas
       filepath: /home/arthur/dev-perso/scikit-learn/sklearn-env/lib/python3.12/site-packages/numpy.libs/libscipy_openblas64_-8fb3d286.so
        version: 0.3.30
threading_layer: pthreads
   architecture: Haswell

       user_api: blas
   internal_api: openblas
    num_threads: 16
         prefix: libscipy_openblas
       filepath: /home/arthur/dev-perso/scikit-learn/sklearn-env/lib/python3.12/site-packages/scipy.libs/libscipy_openblas-b75cc656.so
        version: 0.3.29.dev
threading_layer: pthreads
   architecture: Haswell

       user_api: openmp
   internal_api: openmp
    num_threads: 12
         prefix: libgomp
       filepath: /usr/lib/x86_64-linux-gnu/libgomp.so.1.0.0
        version: None
```

COMMENT:
You're right in that `_weighted_percentile` is not scale invariant for float constants: https://github.com/scikit-learn/scikit-learn/pull/31775#issuecomment-3288112176

For item 1, with average=False (default):

```python
x = np.arange(10)
w = np.ones_like(x)
print(
    _weighted_percentile(x, w),
    _weighted_percentile(x, w * 3 * np.pi),
    _weighted_percentile(x, w * 30 * np.pi)
)
```
 you get: 4.0 4.0 5.0

So it's inaccurate to say it is a problem specifically with `average=True`.

Note that this behaviour is consistent with what would happen with numpy quantile "inverted_cdf" (the same as `average=False`), which is the only method in numpy that supports weights:

```python
x = np.arange(10)
w = np.ones_like(x)
print(
    np.quantile(x, 0.5, method="inverted_cdf", weights=w),
    np.quantile(x, 0.5, method="inverted_cdf", weights=w * 3 * np.pi),
    np.quantile(x, 0.5, method="inverted_cdf", weights=w * 30 * np.pi)
 )
```
gives: 4 4 5

For item 2, yes it is not scale invariant for floats.

For item 3, "average_inverted_cdf" (`average=True`) (see [Hyndman and Fan](https://www.academia.edu/download/51448589/Sample_Quantiles_in_Statistical_Packages20170120-21949-1uaj3m0.pdf)) sets $`m=0`$, such that $`j\le pn\lt j+1`$ (p=percentile, n=sample size), and:

```math
\gamma = \left\{
\begin{array}{ll}
\frac{1}{2} & \text{if } g = 0 \\
1 & \text{if } g > 0
\end{array}
\right.
```
in the general equation $`Q_{i}(p)=(1-\gamma)X_{(j)}+\gamma X_{(j+1)}`$

I suspect what happens is that `g` is never 0 (for large sample size and float weights, its just unlikely for any of cumulative sum of weight values to match exactly 0.5*total cumulative sum) so you always use the 'higher' data point (i.e. $`\gamma`$ is never 1/2)

> So maybe one way to fix this would be to transform/scale the weights in a way that summing them causes no rounding errors. [...]. I will try to explore this idea.

Open to fixes if you have ideas (though with your example of pi, an irrational number, that seems difficult)

> This shows that basically, for large arrays, the option average=True and average=False are basically equivalent.

Honestly, yes. For large arrays I think any quantile method ([numpy](https://numpy.org/doc/stable/reference/generated/numpy.quantile.html) implements 13, incl. 4 'old' ones for backwards compatibility) would essentially be the same, BUT methods do differ in properties, specifically symmetry. "average_inverted_cdf" is the only discrete method that is symmetrical. Most of the continuous (interpolated) methods are symmetrical.

If you're interested there is some discussion here about weights and quantiles: https://github.com/data-apis/array-api-extra/issues/340 and the PR adding weights in numpy https://github.com/numpy/numpy/pull/24254

Labels: Documentation