Repository: scipy/scipy
Pull request: BUG: stats: use mean behavior for percentileofscore in bootstrap
#### Reference issue
Closes gh-15883.

#### What does this implement/fix?
This updates the behavior of the percentileofscore implementation used by bootstrap to use the 'mean' behavior for calculating percentile, which avoids the 'strict' behavior of assigning the min value the zeroth percentile.


Problem statement:
BUG: stats.bootstrap bca implementation triggers ValueError for `min` statistic
### Describe your issue.

`stats.bootstrap` uses a custom implementation of `percentileofscore` that exclusively uses the `strict` behavior from the primary implementation.

As a result, the `min` statistic is always identified as the zeroth percentile, for all datasets. The bca then chooses -inf for z0_hat, which causes a ValueError when the lower percentile is selected as -inf.

I recognize that using bootstrap to estimate confidence intervals of `min` and `max` is not really intended behavior. However, I believe the behavior should be consistent for both of those estimators, and that the `mean` behavior of `percentileofscore` should be used instead.

### Reproducing Code Example

```python
import scipy.stats as stats
import numpy as np
norm = stats.norm.rvs(0, 1, 100)
stats.bootstrap((norm,), statistic=np.min)
```


### Error message

```shell
ValueError                                Traceback (most recent call last)
Input In [1], in <cell line: 4>()
      2 import numpy as np
      3 norm = stats.norm.rvs(0, 1, 100)
----> 4 stats.bootstrap((norm,), statistic=np.min)

File ~/miniconda3/envs/analytics/lib/python3.9/site-packages/scipy/stats/_bootstrap.py:454, in bootstrap(data, statistic, vectorized, paired, axis, confidence_level, n_resamples, batch, method, random_state)
    451         return np.percentile(a=a, q=q, axis=-1)
    453 # Calculate confidence interval of statistic
--> 454 ci_l = percentile_fun(theta_hat_b, interval[0]*100)
    455 ci_u = percentile_fun(theta_hat_b, interval[1]*100)
    456 if method == 'basic':  # see [3]

File ~/miniconda3/envs/analytics/lib/python3.9/site-packages/scipy/stats/_bootstrap.py:78, in _percentile_along_axis(theta_hat_b, alpha)
     76 for indices, alpha_i in np.ndenumerate(alpha):
     77     theta_hat_b_i = theta_hat_b[indices]
---> 78     percentiles[indices] = np.percentile(theta_hat_b_i, alpha_i)
     79 return percentiles[()]

File <__array_function__ internals>:5, in percentile(*args, **kwargs)

File ~/miniconda3/envs/analytics/lib/python3.9/site-packages/numpy/lib/function_base.py:3866, in percentile(a, q, axis, out, overwrite_input, interpolation, keepdims)
   3864 q = asanyarray(q)  # undo any decay that the ufunc performed (see gh-13105)
   3865 if not _quantile_is_valid(q):
-> 3866     raise ValueError("Percentiles must be in the range [0, 100]")
   3867 return _quantile_unchecked(
   3868     a, q, axis, out, overwrite_input, interpolation, keepdims)

ValueError: Percentiles must be in the range [0, 100]
```


### SciPy/NumPy/Python version information

1.7.3 1.21.2 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)

Hints from issue discussion:
I'm not able to reproduce this `ValueError` issue in `main`.  This is the output I get for the provided example:
```
C:\Users\matth\scipy\scipy\stats\_resampling.py:134: RuntimeWarning: invalid value encountered in double_scalars
  alpha_1 = ndtr(z0_hat + num1/(1 - a_hat*num1))
C:\Users\matth\scipy\scipy\stats\_resampling.py:136: RuntimeWarning: invalid value encountered in double_scalars
  alpha_2 = ndtr(z0_hat + num2/(1 - a_hat*num2))
C:\Users\matth\scipy\scipy\stats\_resampling.py:102: BootstrapDegenerateDistributionWarning: The bootstrap distribution is degenerate; the confidence interval is not defined.
  warnings.warn(BootstrapDegenerateDistributionWarning())
Out[74]: BootstrapResult(confidence_interval=ConfidenceInterval(low=nan, high=nan), standard_error=0.11787258252011637)
```
Is there a seed I should use to get the result you got? It's also possible that the particular error here is already fixed (e.g. by gh-14351).

It sounds like that still does not resolve your issue? You would prefer that `percentileofscore` effectively use a different `kind` option, and I'll have to look into that more closely. Are you able to provide more information about why that would be more correct in general, or do you think it should be a user-configurable option? Assuming the `ValueError` is really resolved, would you still consider this to be a _defect_ (as in the algorithm is not implemented correctly), or is this more of an enhancement?
Indeed -- it appears that #14351 prevents the ValueError by directly checking for the `nan` case. 

That said, in this case I do not think you have a degenerate distribution, and I think the confidence interval should be defined, with the `low` value equal to the sample min. 

The existing `scipy.stats` docs for `percentileofscore` describe the `mean` kind as the one often used in testing, which I think is why that would be more correct in general. Further, as a thought experiment, we should not expect the confidence interval calculation to be impossible when we're checking for the `min` statistic but not for the `max` statistic (which is the case with the current `strict` implementation). So I think I would call this a defect (albeit a minor one) rather than an enhancement.

I could also see the case for `rank` method being the appropriate choice, and that has the benefit of being the default of the current `percentileofscore` implementation. I'm not sure how this change would affect existing tests, but would love to create a simple PR if you agree this is a worthy change.
> but would love to create a simple PR if you agree this is a worthy change

Sure. I'd be curious to see what changes in existing tests; hopefully it's not too disruptive for typical cases. I'm sure I originally used `stats.percentileofscore` in the original implementation; I only would have written the custom implementation to improve performance for vectorized input. Could you check how it impacts performance in that case? (For instance, test the example from the documentation "Rather than writing a loop, we can also determine the confidence intervals for all 1000 samples at once.") Probably inconsequential in the scheme of things, but good to check.

It's interesting that you brought this up today. Just this morning I was thinking about gh-15869 and whether something similar would be needed for `bootstrap`. Originally, I thought we would, then this morning I changed my mind because I only remembered how `percentile` was used, forgetting about `percentileofscore`. Now I'm thinking that simply using `percentileofscore` with a different `kind` argument may not be quite enough; you might also need to include a tolerance factor to avoid missing numerically distinct (but theoretically equal) values in the bootstrap distribution.
I'll update here with performance notes, but this had no impact on the tests locally, which I found reassuring.

I'm not sure on the tolerance question--my initial thought is that, even with percentileofscore, we are already allowing for that by the simultaneous `<` and `<=` checks, but I have not read https://github.com/scipy/scipy/pull/15892 closely.

Commit message: Update scipy/stats/tests/test_resampling.py

Patch to review:
diff --git a/scipy/stats/_resampling.py b/scipy/stats/_resampling.py
index 55a99d4edc4e..9bdabec6e46c 100644
--- a/scipy/stats/_resampling.py
+++ b/scipy/stats/_resampling.py
@@ -80,12 +80,13 @@ def _bootstrap_resample(sample, n_resamples=None, random_state=None):
 
 def _percentile_of_score(a, score, axis):
     """Vectorized, simplified `scipy.stats.percentileofscore`.
+    Uses logic of the 'mean' value of percentileofscore's kind parameter.
 
     Unlike `stats.percentileofscore`, the percentile returned is a fraction
     in [0, 1].
     """
     B = a.shape[axis]
-    return (a < score).sum(axis=axis) / B
+    return ((a < score).sum(axis=axis) + (a <= score).sum(axis=axis)) / (2 * B)
 
 
 def _percentile_along_axis(theta_hat_b, alpha):
diff --git a/scipy/stats/tests/test_resampling.py b/scipy/stats/tests/test_resampling.py
index 9b733d774719..e8853fd2267e 100644
--- a/scipy/stats/tests/test_resampling.py
+++ b/scipy/stats/tests/test_resampling.py
@@ -388,6 +388,23 @@ def test_bootstrap_gh15678(method):
     assert isinstance(res.standard_error, np.float64)
 
 
+def test_bootstrap_min():
+    # Check that gh-15883 is fixed: percentileofscore should
+    # behave according to the 'mean' behavior and not trigger nan for BCa
+    rng = np.random.default_rng(1891289180021102)
+    dist = stats.norm(loc=2, scale=4)
+    data = dist.rvs(size=100, random_state=rng)
+    true_min = np.min(data)
+    data = (data,)
+    res = bootstrap(data, np.min, method="BCa", n_resamples=100,
+                    random_state=np.random.default_rng(3942))
+    assert true_min == res.confidence_interval.low
+    res2 = bootstrap(-np.array(data), np.max, method="BCa", n_resamples=100,
+                     random_state=np.random.default_rng(3942))
+    assert_allclose(-res.confidence_interval.low, res2.confidence_interval.high)
+    assert_allclose(-res.confidence_interval.high, res2.confidence_interval.low)
+
+
 def test_jackknife_resample():
     shape = 3, 4, 5, 6
     np.random.seed(0)


Code review comments:
diff_hunk: @@ -388,6 +388,23 @@ def test_bootstrap_gh15678(method):
     assert isinstance(res.standard_error, np.float64)
 
 
+def test_bootstrap_min():
+    # Check that gh-15883 is fixed: percentileofscore should
+    # behave according to the 'mean' behavior and not trigger nan for BCa
+    rng = np.random.default_rng(1891289180021102)
+    dist = stats.norm(loc=2, scale=4)
+    data = dist.rvs(size=100, random_state=rng)
+    true_min = np.min(data)
+    data = (data,)
+    res = bootstrap(data, np.min, method="BCa", n_resamples=100,
+                    random_state=np.random.default_rng(3942))
+    assert true_min == res.confidence_interval.low
+    res2 = bootstrap(-np.array(data), np.max, method="BCa", n_resamples=100,
+                     random_state=np.random.default_rng(3942))
+    assert_allclose(-res.confidence_interval.low, res2.confidence_interval.high)
+    assert_allclose(-res.confidence_interval.high, res2.confidence_interval.low)
path: scipy/stats/tests/test_resampling.py
text: @user1:
```suggestion
    assert_allclose(-res.confidence_interval.low,
                    res2.confidence_interval.high)
    assert_allclose(-res.confidence_interval.high,
                    res2.confidence_interval.low)
```