Repository: scipy/scipy
Pull request: ENH: stats.bootstrap: add BCa method for multi-sample statistic
#### Reference issue
Closes gh-16433
Addresses item in gh-13371

#### What does this implement/fix?
Currently, `stats.bootstrap` only supports `method='BCa'` for single-sample statistics. This adds support for multi-sample statistics.

I need to add a test. @FlorinAndrei can you help me devise one? One way to do this would be to provide output from R's `boot.ci`  to compare against. Ideally, this would be performed with many (say ~100k resamples) to avoid discrepancies due to the choice of random samples, and it would be nice if the data/statistic chosen result in a BCa interval that is substantially different from the `'perc'` and `'basic'` versions. Also, if you could provide the bootstrap replicates of the statistic from R (`boot.out$t`), I could check the BCa interval that I'd get with the same bootstrap distribution. There might also be other functions in the `boot` package that could provide relevant output to compare against.

#### Additional information
We would have supported multi-sample statistics in the original PR, but we were unsure how the jackknife should be performed (to approximate the acceleration parameter) when there are multiple samples.

~~In this PR, the jackknife removes one observation from _one_ sample at a time because this seems to be the recommendation in [Missing Data, Imputation, and the Bootstrap](https://www.jstor.org/stable/2290846?seq=1). This is in contrast to what would seem like the more careful (and _much_ more time-consuming) option: removing one observation from _every_ sample (in each possible way), which is the recommendation of [Quantifying effects in two-sample environmental experiments using bootstrap confidence intervals](https://www.sciencedirect.com/science/article/abs/pii/S1364815205002239?casa_token=OcSrOn8nclMAAAAA:FKgCGSdCysMAjQl6Qf7v58encmDQmdSgOdV3iv1uFI48ePUgbz7HjpzaYBKMF9bRs7v4Ll0OGcc). However, that article mentions performing a comparison with the fast method implemented here, and they found "only very small differences (vs. nested sums) in case of robust measures and small differences in case of non-robust measures, not affecting the main conclusions of this paper". I'll continue to look into the theory, but @steppi if you know or are interested in the theory of the bootstrap, this is (yet) another one I could use some help with!~~

Efron presents the calculation of the acceleration parameter in the multi-sample case in Chapter 15.4 of An Introduction to the Bootstrap, specifically Equation 15.36.

See https://github.com/scipy/scipy/issues/16433#issuecomment-1161358545 for additional discussion.

Problem statement:
ENH: scipy.stats.bootstrap() should do BCa for multivariate statistics too
### Is your feature request related to a problem? Please describe.

I've seen this several times while looking at datasets. In R I can do BCa correction any time. When I try to do the same analysis with Python and I want to apply the BCa correction, it only works for univariate.

### Describe the solution you'd like.

Would be really nice to have BCa available for all analyses.

### Describe alternatives you've considered.

R-lang - I mean, it works already.

### Additional context (e.g. screenshots, GIFs)

_No response_

Hints from issue discussion:
Ok, i can add it. The reason we didn't is that there is not much literature that we could find about a multi-sample jackknife. It is conceptually simple to generalize, but the algorithm is time-consuming. So we started thinking about using bootstrap within bootstrap... and decided to let it go for the time being.
Without looking at the R code, do you have a reference for doing the jackknife part with multiple samples? (We found one, but it didn't seem like a super reliable reference.)
I don't have it, sorry. All I know is that R does BCa just fine in either case. No idea how it does that, it's probably in the code.
multivariate or multi-sample ?

those are different things (in most applications)

> multivariate or multi-sample ?

Let me provide a concrete example to make it clear what I want to do.

[BirdPecks.csv](https://github.com/scipy/scipy/files/8943744/BirdPecks.csv)

This code fails:

```python
import pandas as pd
import numpy as np
from scipy import stats

BirdPecks = pd.read_csv("BirdPecks.csv")

pg1 = BirdPecks[BirdPecks["group"] == 1]["pecks"].to_numpy()
pg2 = BirdPecks[BirdPecks["group"] == 2]["pecks"].to_numpy()

bdist = []
global bdist
Nboot = 10000

def boot_mean(x, y):
    global bdist
    m = np.mean(x) - np.mean(y)
    bdist.append(m)
    return m

bs_out = stats.bootstrap(
    data=(pg1, pg2),
    statistic=boot_mean,
    vectorized=False,
    paired=False,
    confidence_level=0.9,
    n_resamples=Nboot,
    random_state=123,
    method="BCa"
)
```
Yup. Looks like you'd also like the bootstrap distribution to be returned, based on `bdist`? What would you like to use it for?
Here is why it wasn't returned originally: https://github.com/scipy/scipy/pull/13371#issuecomment-833176091
Suggest how you'd use it, and I'd be happy to return it.

And here's the original discussion about multi-sample statistics: https://github.com/scipy/scipy/pull/13371#issuecomment-809151181. 
What R function do you use to do this, BTW? (e.g. If you use `boot`, do you use `strata` to distinguish between the two samples?) 

I looked at the [documentation](https://cran.r-project.org/web/packages/bcaboot/bcaboot.pdf) and corresponding [paper ](https://efron.ckirby.su.domains/papers/2018Automatic-Construction-BCIs.pdf) of Efron and Narasimhan's package `bcaboot`. There are some cool ideas there for extending `bootstrap`, but multi-sample statistics are explicitly not supported:

![image](https://user-images.githubusercontent.com/6570539/174729835-f00ae970-c41d-45e6-82cd-90193fa683dc.png)

For future reference to others who find this issue: I believe that _multivariate_ statistics are already handled: statistics of multivariate samples (i.e. multi-input statistics) are supported with `paired=True`, and vector-valued (multi-output) statistics are supported as of gh-15491. (I can imagine ways in which all this could be generalized - e.g. multiple multivariate samples - but, but let's take it one step at a time.)
> Looks like you'd also like the bootstrap distribution to be returned, based on bdist? What would you like to use it for?

That's not always needed, of course, but you could run tests on the bootstrap distribution, compare it with known distributions (to what extent does it diverge from those, test a hypothesis, etc), assess central values, etc. It greatly extends the usefulness of the original function, besides the basic test you do with the `statistic` parameter. There are times when just graphing it is useful.

It could be an optional parameter, return the distribution only when needed.

> What R function do you use to do this, BTW?

Data set:

[repair.csv](https://github.com/scipy/scipy/files/8953607/repair.csv)

R code:

```r
rp <- read.csv('repair.csv')
rp$type <- factor(rp$type)

require(boot)
bootMeanDiff <- function(d, i){
  means <- tapply(d[i, 1], d[, 2], mean)
  means[1] - means[2]
}

boot.object <- boot(rp, bootMeanDiff, R = 5000, strata = rp$type)
bca.ci = boot.ci(boot.object, conf=.95, type='bca')$bca[4:5]
print(bca.ci)
```
There are also libraries that buid on top of Scipy, and having the bootstrap distribution returned directly from the main Scipy function call would simplify things a lot (as opposed to doing tricks with `bd.append()` like I did).

Commit message: ENH: stats.bootstrap: add BCa method for multi-sample statistic

Patch to review:
diff --git a/scipy/stats/_resampling.py b/scipy/stats/_resampling.py
index d1c938f28d4f..6e5e2d867e1a 100644
--- a/scipy/stats/_resampling.py
+++ b/scipy/stats/_resampling.py
@@ -103,17 +103,27 @@ def _percentile_along_axis(theta_hat_b, alpha):
 def _bca_interval(data, statistic, axis, alpha, theta_hat_b, batch):
     """Bias-corrected and accelerated interval."""
     # closely follows [2] "BCa Bootstrap CIs"
-    sample = data[0]  # only works with 1 sample statistics right now
+    # for extension to multiple samples, see:
+    # https://github.com/scipy/scipy/issues/16433#issuecomment-1161358545
 
     # calculate z0_hat
-    theta_hat = np.asarray(statistic(sample, axis=axis))[..., None]
+    theta_hat = np.asarray(statistic(*data, axis=axis))[..., None]
     percentile = _percentile_of_score(theta_hat_b, theta_hat, axis=-1)
     z0_hat = ndtri(percentile)
 
     # calculate a_hat
     theta_hat_i = []  # would be better to fill pre-allocated array
-    for jackknife_sample in _jackknife_resample(sample, batch):
-        theta_hat_i.append(statistic(jackknife_sample, axis=-1))
+    for i, sample in enumerate(data):
+        # jackknife resampling will add an axis prior to the last axis that
+        # corresponds with the different jackknife resamples. Do the same for
+        # each sample of the data to ensure broadcastability. We need to
+        # create a copy of the list containing the samples anyway, so do this
+        # in the loop to simplify the code. This is not the bottleneck...
+        samples = [np.expand_dims(sample, -2) for sample in data]
+        for jackknife_sample in _jackknife_resample(sample, batch):
+            samples[i] = jackknife_sample
+            samples = _broadcast_arrays(samples, axis=-1)
+            theta_hat_i.append(statistic(*samples, axis=-1))
     theta_hat_i = np.concatenate(theta_hat_i, axis=-1)
     theta_hat_dot = theta_hat_i.mean(axis=-1, keepdims=True)
     num = ((theta_hat_dot - theta_hat_i)**3).sum(axis=-1)
@@ -199,10 +209,6 @@ def statistic(i, axis=-1, data=data_iv, unpaired_statistic=statistic):
     if method not in methods:
         raise ValueError(f"`method` must be in {methods}")
 
-    message = "`method = 'BCa' is only available for one-sample statistics"
-    if not paired and n_samples > 1 and method == 'bca':
-        raise ValueError(message)
-
     random_state = check_random_state(random_state)
 
     return (data_iv, statistic, vectorized, paired, axis_int,
@@ -279,10 +285,9 @@ def bootstrap(data, statistic, *, n_resamples=9999, batch=None,
         The confidence level of the confidence interval.
     method : {'percentile', 'basic', 'bca'}, default: ``'BCa'``
         Whether to return the 'percentile' bootstrap confidence interval
-        (``'percentile'``), the 'reverse' or the bias-corrected and accelerated
-        bootstrap confidence interval (``'BCa'``).
-        Note that only ``'percentile'`` and ``'basic'`` support multi-sample
-        statistics at this time.
+        (``'percentile'``), the 'basic' (AKA 'reverse') bootstrap confidence
+        interval (``'basic'``), or the bias-corrected and accelerated bootstrap
+        confidence interval (``'BCa'``).
     random_state : {None, int, `numpy.random.Generator`,
                     `numpy.random.RandomState`}, optional
 
diff --git a/scipy/stats/tests/test_resampling.py b/scipy/stats/tests/test_resampling.py
index 2e3d053f4405..36a790059cc1 100644
--- a/scipy/stats/tests/test_resampling.py
+++ b/scipy/stats/tests/test_resampling.py
@@ -59,16 +59,6 @@ def test_bootstrap_iv():
     with pytest.raises(ValueError, match=message):
         bootstrap(([1, 2, 3],), np.mean, method='ekki')
 
-    message = "`method = 'BCa' is only available for one-sample statistics"
-
-    def statistic(x, y, axis):
-        mean1 = np.mean(x, axis)
-        mean2 = np.mean(y, axis)
-        return mean1 - mean2
-
-    with pytest.raises(ValueError, match=message):
-        bootstrap(([.1, .2, .3], [.1, .2, .3]), statistic, method='BCa')
-
     message = "'herring' cannot be used to seed a"
     with pytest.raises(ValueError, match=message):
         bootstrap(([1, 2, 3],), np.mean, random_state='herring')


Code review comments:
diff_hunk: @@ -103,17 +103,27 @@ def _percentile_along_axis(theta_hat_b, alpha):
 def _bca_interval(data, statistic, axis, alpha, theta_hat_b, batch):
     """Bias-corrected and accelerated interval."""
     # closely follows [2] "BCa Bootstrap CIs"
-    sample = data[0]  # only works with 1 sample statistics right now
+    # for extension to multiple samples, see:
+    # https://github.com/scipy/scipy/issues/16433#issuecomment-1161358545
 
     # calculate z0_hat
-    theta_hat = np.asarray(statistic(sample, axis=axis))[..., None]
+    theta_hat = np.asarray(statistic(*data, axis=axis))[..., None]
     percentile = _percentile_of_score(theta_hat_b, theta_hat, axis=-1)
     z0_hat = ndtri(percentile)
 
     # calculate a_hat
     theta_hat_i = []  # would be better to fill pre-allocated array
-    for jackknife_sample in _jackknife_resample(sample, batch):
-        theta_hat_i.append(statistic(jackknife_sample, axis=-1))
+    for i, sample in enumerate(data):
+        # jackknife resampling will add an axis prior to the last axis that
+        # corresponds with the different jackknife resamples. Do the same for
+        # each sample of the data to ensure broadcastability. We need to
+        # create a copy of the list containing the samples anyway, so do this
+        # in the loop to simplify the code. This is not the bottleneck...
+        samples = [np.expand_dims(sample, -2) for sample in data]
+        for jackknife_sample in _jackknife_resample(sample, batch):
+            samples[i] = jackknife_sample
+            samples = _broadcast_arrays(samples, axis=-1)
+            theta_hat_i.append(statistic(*samples, axis=-1))
     theta_hat_i = np.concatenate(theta_hat_i, axis=-1)
     theta_hat_dot = theta_hat_i.mean(axis=-1, keepdims=True)
     num = ((theta_hat_dot - theta_hat_i)**3).sum(axis=-1)
path: scipy/stats/_resampling.py
text: @author:
These are not quite right. The individual terms from different samples should not be weighted equally.

diff_hunk: @@ -103,17 +103,27 @@ def _percentile_along_axis(theta_hat_b, alpha):
 def _bca_interval(data, statistic, axis, alpha, theta_hat_b, batch):
     """Bias-corrected and accelerated interval."""
     # closely follows [2] "BCa Bootstrap CIs"
-    sample = data[0]  # only works with 1 sample statistics right now
+    # for extension to multiple samples, see:
+    # https://github.com/scipy/scipy/issues/16433#issuecomment-1161358545
 
     # calculate z0_hat
-    theta_hat = np.asarray(statistic(sample, axis=axis))[..., None]
+    theta_hat = np.asarray(statistic(*data, axis=axis))[..., None]
     percentile = _percentile_of_score(theta_hat_b, theta_hat, axis=-1)
     z0_hat = ndtri(percentile)
 
     # calculate a_hat
     theta_hat_i = []  # would be better to fill pre-allocated array
-    for jackknife_sample in _jackknife_resample(sample, batch):
-        theta_hat_i.append(statistic(jackknife_sample, axis=-1))
+    for i, sample in enumerate(data):
+        # jackknife resampling will add an axis prior to the last axis that
+        # corresponds with the different jackknife resamples. Do the same for
+        # each sample of the data to ensure broadcastability. We need to
+        # create a copy of the list containing the samples anyway, so do this
+        # in the loop to simplify the code. This is not the bottleneck...
+        samples = [np.expand_dims(sample, -2) for sample in data]
+        for jackknife_sample in _jackknife_resample(sample, batch):
+            samples[i] = jackknife_sample
+            samples = _broadcast_arrays(samples, axis=-1)
path: scipy/stats/_resampling.py
text: @author:
`samples` needs a different name in case the last batch is a different length than the rest.