Repository: scipy/scipy
Pull request: ENH: stats: add general permutation hypothesis test
#### Reference issue
<!--Example: Closes gh-WXYZ.-->
gh-4824 - this PR is a generalization of that PR to a user-provided statistic of any number of samples
gh-4399 - probably closes this (see "Additional Information")
gh-4519 - supersedes that PR (see "Additional Information")
gh-13263 - in a followup PR, this function should be used to replace the exact statistic implementation. (A `permutation` argument could be added.) For now, that exact statistic calculation is used as a test.
gh-13371 - another resampling method (bootstrapping)
~~gh-13661 - this PR should wait for that PR to merge so they can share some code.~~
~~gh-14206 - this PR could wait for that PR to merge to reduce some code duplication~~

#### What does this implement/fix?
<!--Please explain your changes.-->
This adds `permutation_test`, a function for performing a permutation test for user-provided statistics. It performs a randomized permutation test when the specified number of `permutations` is less than number of possible "permutations" of the data, otherwise it performs an exact test. 

Three types of permutations can be performed:
* `'samples'` - observations are permuted between samples without changing their pairing with observations in other samples. This type of permutation is consistent with the null hypothesis of tests like the paired-sample t-test and the Wilcoxon signed rank test.
* `'pairings'` - observations are permuted between pairings without changing the sample to which they are assigned. This type of permutation is consistent with the null hypothesis of many association tests (e.g. Pearson, Spearman, Kendall, Page).
* `'both'` - both types of permutation are performed. This is consistent with the null hypothesis of tests like the independent sample t-test, the Mann-Whitney U test, one-way ANOVA, and the Kruskal-Wallis H test.

With appropriate definition of the test statistic and choice of pairing, this function can reproduce the results of most of the statistical tests in SciPy.

#### Benefits
- Exact tests for small samples, even when there are ties
- More accurate tests for moderate samples when the asymptotic approximation is not appropriate and an exact test is not implemented
- Easily perform hyothesis tests not yet implemented in SciPy
- Conveniently check the accuracy of (new) exact hypothesis test implementations

#### Examples
Perform an `exact` Mann-Whitney U test using `permutation_test` and compare to exact p-value produced by `mannwhitneyu`:

```python3
import numpy as np
from scipy import stats
from scipy.stats import permutation_test

np.random.seed(0)
alternative = 'less'
x = stats.norm.rvs(size=8)
y = stats.norm.rvs(size=9, scale=0.75)

expected = stats.mannwhitneyu(x, y, alternative=alternative, method='exact')

def statistic(x, y, axis):
    return stats.mannwhitneyu(x, y, alternative=alternative, axis=axis).statistic

res = permutation_test((x, y), statistic, alternative=alternative, vectorized=True)

print(f'Permutation p-value: {res.pvalue}')
print(f'Expected p-value: {expected.pvalue}')
print(f'Difference: {res.pvalue - expected.pvalue}')
# Permutation p-value: 0.8617030028794734
# Expected p-value: 0.8617030028794733
# Difference: 1.1102230246251565e-16
```

Compare the exact Ansari-Bradley null distribution for the given tie pattern against the asymptotic approximation:
<details>

```python3
import matplotlib.pyplot as plt

x = np.array([4, 10, 3, 1, 7, 2, 6, 7])
y = np.array([5, 8, 6, 2, 10, 3, 9, 4])

def statistic(x, y):
    return stats.ansari(x, y).statistic

res = permutation_test((x, y), statistic)  # considers ties
bins = np.linspace(np.min(res.null_distribution), np.max(res.null_distribution), 50)
pdf = stats.norm.pdf(bins, loc=36, scale=np.sqrt(21.6))  # from inside `ansari`

plt.hist(res.null_distribution, bins=bins, density=True)
plt.plot(bins, pdf)
plt.legend(("Exact Null Distribution, Considering Ties",
            "Asymptotic Approximation"))
```

</details>

![image](https://user-images.githubusercontent.com/6570539/115627456-df0dc780-a2b3-11eb-838f-3e20754427da.png)


Use `permutation_test` to perform a correlation test using 10000 permutations, and compare against asymptotic p-value from `spearmanr`.

<details>

```python3
x = stats.norm.rvs(size=500)
y = stats.norm.rvs(size=500)

alternative = 'greater'
expected = stats.spearmanr(x, y, alternative=alternative)

def statistic(x):
    return stats.spearmanr(x, y).correlation

res = permutation_test((x,), statistic, permutation_type='pairings',
                       permutations=10000, alternative=alternative)

print(f'Relative Difference: '
      f'{np.abs(res.pvalue-expected.pvalue)/expected.pvalue}')
# Relative Difference: 0.00010797224590931033
```

</details>

Use `permutation_test` to perform an exact one-sample Wilcoxon signed-rank , and compare against exact p-value from `wilcoxon`.
<details>

```python3
import numpy as np
from scipy import stats
from scipy.stats import permutation_test

np.random.seed(0)
alternative = 'less'
x = stats.norm.rvs(size=8)
y = stats.norm.rvs(size=8, scale=0.75)

expected = stats.wilcoxon(x - y, alternative=alternative, mode='exact')

def statistic(z):
    # always use 'less' statistic, otherwise 2nd statistic could be returned
    return stats.wilcoxon(z, alternative='less').statistic

res = permutation_test((x - y,), statistic, permutation_type='samples',
                       alternative=alternative)

print(f'Permutation p-value: {res.pvalue}')
print(f'Expected p-value: {expected.pvalue}')
print(f'Difference: {res.pvalue - expected.pvalue}')
# Permutation p-value: 0.07421875
# Expected p-value: 0.07421875
# Difference: 0.0
```

</details>

Use `permutation_test` to perform an exact three-sample Page trend test , and compare against exact p-value from `page_trend_test`.

<details>

```python3
import numpy as np
from scipy import stats
from scipy.stats import permutation_test

np.random.seed(0)

x = stats.norm.rvs(size=4)
y = stats.norm.rvs(size=4)
z = stats.norm.rvs(size=4)

expected = stats.page_trend_test(np.vstack([x, y, z]), method='exact')

def statistic(*data):
    # always use 'less' statistic, otherwise 2nd statistic could be returned
    return stats.page_trend_test(data).statistic

res = permutation_test((x, y, z,), statistic, permutation_type='pairings',
                       alternative='greater')

print(f'Permutation p-value: {res.pvalue}')
print(f'Expected p-value: {expected.pvalue}')
print(f'Difference: {res.pvalue - expected.pvalue}')
# Permutation p-value: 0.07421875
# Expected p-value: 0.07421875
# Difference: 0.0
```

</details>

#### Additional information
<details>

<!--Any additional information you think is important.-->
In gh-13661, @jb-leger and I have been discussing enhancements that will improve the speed.
~~We can generalize this to n-samples if desired. We should decide before merging whether we'll want to make that part of this function, because we'd either want to take in all `data` as a tuple, rather than separate `x` and `y` arguments, or use `*args`.~~ _Update: done; `data` must be a sequence._
~~We should discuss the definition of two-sided p-values that should be used.~~ _Done. Here we use twice the one-sided p-value. This is a different convention from permutation t-test. As an enhancement, this function could support several conventions for two-sided p-values._
~~This PR should wait for gh-13661 to merge first. We should combine a lot of the code.~~ _Done. Also reused some code from bootstrap._
~~This should also wait for gh-4933 to merge. We can test this function against the MWU null distribution.~~ _Done._
</details>

Regarding gh-4399 and gh-4519:
I compared the performance of this PR against the code in gh-4519 for two-sample ANOVA, as this PR had not yet been generalized to n-samples. 

```
>>> def my_statistic(*data, axis):
...     return stats.f_oneway(*data, axis=axis)[0]
>>> np.random.seed(1)
>>> x = np.random.rand(100)
>>> y = np.random.rand(100)
>>> %timeit f_oneway(x, y, permutations=10000)  # gh-4519
684 ms ± 5.32 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>>> %timeit stats.permutation_test((x, y), my_statistic, vectorized=True, permutations=10000, alternative='greater')  # this PR
90.5 ms ± 929 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

<details>

For larger samples:
```
>>> np.random.seed(1)
>>>  x = np.random.rand(1000)
>>>  y = np.random.rand(1000)
>>>  %timeit f_oneway(x, y, permutations=10000)
1.91 s ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>>>  %timeit stats.permutation_test((x, y), my_statistic, vectorized=True, permutations=10000, alternative='greater')
692 ms ± 16.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

For a smaller sample with more permutations:
```
>>> np.random.seed(1)
>>>  x = np.random.rand(10)
>>>  y = np.random.rand(10)
>>> %timeit stats.permutation_test((x, y), my_statistic, vectorized=True, permutations=np.inf, alternative='greater')  # exact test
1.35 s ± 10.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>>> %timeit f_oneway(x, y, permutations=184756)  # same number of random permutations (not exact)
10.5 s ± 31.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

</details>

As this PR is generalized to user-provided statistics and N-d array samples (not just 2D), I'd suggest that this PR closes gh-4399. Now that it has been extended to multi-sample statistics, gh-4519 can also be closed.

I also compared against the two-sample example of the [`permute` package](https://statlab.github.io/permute/user/two-sample.html).
```python3
%timeit two_sample(maleid, femaleid, stat='t', alternative='two-sided', seed=20, reps=10000)  # permute package
# 3.42 s ± 10.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```
vs
```python3
from scipy.stats import ttest_ind
maleid = [4, 4, 5, 5, 5, 4, 4, 5, 4, 3, 5, 4]
femaleid = [1, 1, 5, 4, 4, 5, 4, 5]
def statistic(x, y, axis):
    return ttest_ind(x, y, axis=axis)[0]
%timeit permutation_test((maleid, femaleid), statistic, vectorized=True, permutations=10000)  # this PR
# 28.2 ms ± 40.5 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```
For that data, this PR can also perform an exact test:
```python3
%timeit permutation_test((maleid, femaleid), statistic, vectorized=True)
# 916 ms ± 5.55 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```
_Update: times increased ~20% after generalizing to n-sample statistics. But let's focus on correctness for now, and we can optimize later. It's already much faster than these alternatives._

Problem statement:
Incorporating in fast permutation tests
Recently, I've been developing fast permutation tests that can make use of numpy/scipy matrix operations.

So far, I have developed 3 permutation tests: a mean test, a t-test and a pearson correlation test.

To get an idea for the timings, I've benchmarked the mean test to a naive implementation on a 60x60 table and  with 1000 permutations. Below are the timings

Naive time [s]: 1.91008838018
Numpy time [s]: 0.0174480279287

If you guys are down, I can start committing these algorithms into scipy. In the near future, I plan on implementing ANOVA tests such as an F-test and a Krustal-Wallis test.

Hints from issue discussion:
@mortonjt I don't quite understand what the difference is between the tests with the same names that are already in `scipy.stats`. Can you elaborate? And do you have them up on Github somewhere already?

I really appreciate the response.

I've add the code a simple test case to the following code
https://github.com/mortonjt/scipy/blob/master/scipy/stats/stats.py
https://github.com/mortonjt/scipy/blob/master/scipy/stats/tests/test_stats.py

The basic idea is it'll calculate a t-test using permutations rather than assuming normality.
In addition, it can calculate multiple t-tests at once for multiple vectors of the same length given a categorical vector of the same length. By making use of numpy matrix multiplication, this is much much faster than looping through and calculating t-statistics for each individual vector.

If the math isn't quite making sense, I'll be happy to add some latex markup in the documentation.

```
def ttest_perm(a, b, axis=0, permutations = 1000):
    """
    Calculates the Welch's T-test for the means of TWO INDEPENDENT samples of scores
    using permutation methods
    This test is an equivalent to scipy.stats.ttest_ind, except it doesn't require the
    normality assumption since it uses a permutation test.
    Parameters
    ----------
    a : array_like
        The array must be 2 dimensional
    b: array_like
        This array must be 1 dimensional and have the same size as
        the length of the specified axis
    axis : int, optional
        Axis can equal None (ravel array first), or an integer (the axis
        over which to operate on a and b).
    permutations: int
        Number of permutations used to calculate p-value
    Returns
    -------
    t : float or array
        The calculated t-statistic.
    prob : float or array
        The two-tailed p-value.
    """
    assert len(a.shape) == 2
    r, c = a.shape
    mat = np.matrix(a)
    if axis == 0: # Test with respect to rows
        assert len(b) == c
    else:
        assert len(b) == r
        mat = mat.transpose()

    perms = _init_categorical_perms(b, permutations=1000)
    perms = np.matrix(perms)
    num_cats = 2
    _, c = perms.shape
    permutations = (c - num_cats) / num_cats

    ## Perform matrix multiplication on data matrix
    ## and calculate sums and squared sums
    _sums  = mat * perms
    _sums2 = np.multiply(mat, mat) * perms

    ## Calculate means and sample variances
    tot =  perms.sum(axis = 0)
    _avgs  = _sums / tot
    _avgs2 = _sums2 / tot
    _vars  = _avgs2 - np.multiply(_avgs, _avgs)
    _samp_vars =  np.multiply(tot, _vars) / (tot-1)

    ## Calculate the t statistic
    idx = np.arange(0, (permutations+1)*num_cats, num_cats, dtype=np.int32)
    denom  = np.sqrt(_samp_vars[:, idx+1] / tot[:,idx+1]  + _samp_vars[:, idx] / tot[:,idx])
    t_stat = np.divide(_avgs[:, idx] - _avgs[:, idx+1], denom)
    ## Calculate the p-values
    cmps =  abs(t_stat[:,1:]) >= abs(t_stat[:,0])
    pvalues = (cmps.sum(axis = 1) + 1.)/(permutations + 1.)

    return map(np.array, map(np.ravel, [t_stat[:, 0],pvalues]))
```

@josef-pkt @jseabold you guys have a better overview of the t-test jungle and whether this fits in `scipy.stats`. What do you think?

If we want to add all the proposed tests, we'll also need to think of a sane API, to prevent the many `ttest_xxx` functions situation for the other ones.

@mortonjt for a 60x60 matrix the calculation time isn't too bad, but how does it scale with rows/cols?

You should avoid using `np.matrix`, using normal arrays (and matrix multiplication with `np.dot`) is very much preferable.

It actually scales pretty well.  

Here are some timings that I have conducted on a 12 core Xeon processor (using multithreaded OpenBLAS)

Here are some of the timings scaling the number of permutations.  N = rows, M = columns

| N | M | Permutations | ttest_ind | numpy |
| --- | --- | --- | --- | --- |
| 500 | 500 | 300 | 13.15 | 0.18 |
| 500 | 500 | 600 | 26.42 | 0.36 |
| 500 | 500 | 900 | 39.21 | 0.52 |

And here are some timings scaling the rows (N).  

| N | M | Permutations | Numpy |
| --- | --- | --- | --- |
| 1000 | 1000 | 5000 | 26.902 |
| 2000 | 1000 | 5000 | 53.137 |
| 3000 | 1000 | 5000 | 79.091 |

Note: I dropped the timings with ttest_ind because they were taking a bit long to complete

This looks really interesting.

> If we want to add all the proposed tests, we'll also need to think of a sane API, to prevent the many ttest_xxx functions situation for the other ones.

couldn't one add a paremter to switch between classical and permutation variants instead of inflating the API?

Alexandre Gramford has written permutation t-test for mne
https://github.com/mne-tools/mne-python/blob/master/mne/stats/permutations.py
that I was targeting for inclusion in statsmodels

Since this is a new area for scipy.stats, it requires some upfront design decisions, especially given the deprecation policy of scipy.  I don't have a strong opinion yet.
- permutation and bootstrap test will be necessary for python either in scipy.stats or statsmodels or both
- efficient permutation needs to be build into the specific functions to make it efficient, e.g. save on unnecessary repeated calculations. This should be considerably faster than permutation or bootstrap "from the outside". I have an example script for another test somewhere in a scipy issue where we had this discussion before. (see motivation of this issue)
- I think the API could avoid inflating the number of functions by adding a `use_permutation` (boolean) and `permutation_kwds` to the current functions for hypothesis tests. But then it would need to offer the other options also. (Actual implementation could also be in a separate private function.)
- Sticking to the current pattern of returning (statistic, pvalue) limits the possibility of providing additional returns which might reuse existing bootstrap or permutation results

for statsmodels
- this will be added -- eventually
- I haven't worked on statsmodels.stats in a while except for hypothesis tests related to estimation models and for maintenance. Unless there is a PR for statsmodels, I would expect this to get added in two versions or so.
- The advantage of statsmodels is that we are not restricted to simple functions with simple returns. Design and refactoring are, however, still unfinished.
- I don't mind code duplication between statsmodels and scipy.stats in some areas. scipy.stats with the basic functions and statsmodels with the "verbose" functions and with extra frills. (I was thinking of copying and duplicating some statistical tests from scipy.stats once we have decided on a "verbose" return/results class for basic hypothesis tests.) 
- statsmodels has joblib as optional dependency which makes it easier to offer a parallel option, same as Alexandre's version.
- extending permutation and bootstrap beyond the simple cases requires supporting infrastructure (for example iterators as for cross-validation in scikit-learn)

In general, t-test is pretty robust to different distributions, some of the other hypothesis test would benefit more from a permutation version. The same pattern should work for all the two or k-sample tests.

To avoid a vaporware issue: Eventually having it in statsmodels doesn't prevent from adding it to scipy.stats. (But I would also be happy to review a PR for statsmodels. :)

correction: The permutation t-test in the above mne link is for one-sample or paired sample. I thought Alexandre and I also discussed two sample tests but maybe I don't remember correctly.

Commit message: DOC: stats: fixes/improvements to docstring per review

[skip azp] [skip actions]

Patch to review:
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index f410f5410d25..cbd02bf331c0 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -301,6 +301,7 @@
    combine_pvalues
    jarque_bera
    page_trend_test
+   permutation_test
 
 .. autosummary::
    :toctree: generated/
diff --git a/scipy/stats/_bootstrap.py b/scipy/stats/_bootstrap.py
index 134ea599fff3..d87f4fa4c532 100644
--- a/scipy/stats/_bootstrap.py
+++ b/scipy/stats/_bootstrap.py
@@ -4,7 +4,7 @@
 from scipy._lib._util import rng_integers
 from dataclasses import make_dataclass
 from ._common import ConfidenceInterval
-from scipy.stats.stats import _broadcast_concatenate
+import scipy.stats.stats as statsstats
 
 
 def _vectorize_statistic(statistic):
@@ -14,7 +14,7 @@ def _vectorize_statistic(statistic):
     def stat_nd(*data, axis=0):
         lengths = [sample.shape[axis] for sample in data]
         split_indices = np.cumsum(lengths)[:-1]
-        z = _broadcast_concatenate(data, axis)
+        z = statsstats._broadcast_concatenate(data, axis)
 
         def stat_1d(z):
             data = np.split(z, split_indices)
@@ -265,15 +265,15 @@ def bootstrap(data, statistic, *, vectorized=True, paired=False, axis=0,
     random_state : {None, int, `numpy.random.Generator`,
                     `numpy.random.RandomState`}, optional
 
-        If `seed` is ``None`` (or `np.random`), the `numpy.random.RandomState`
-        singleton is used.
-        If `seed` is an int, a new ``RandomState`` instance is used,
-        seeded with `seed`.
-        If `seed` is already a ``Generator`` or ``RandomState`` instance then
-        that instance is used.
-
         Pseudorandom number generator state used to generate resamples.
 
+        If `random_state` is ``None`` (or `np.random`), the
+        `numpy.random.RandomState` singleton is used.
+        If `random_state` is an int, a new ``RandomState`` instance is used,
+        seeded with `random_state`.
+        If `random_state` is already a ``Generator`` or ``RandomState``
+        instance then that instance is used.
+
     Returns
     -------
     res : BootstrapResult
diff --git a/scipy/stats/_hypotests.py b/scipy/stats/_hypotests.py
index 897c42ff8ff6..4b6af594f46c 100644
--- a/scipy/stats/_hypotests.py
+++ b/scipy/stats/_hypotests.py
@@ -2,17 +2,20 @@
 from dataclasses import make_dataclass
 import numpy as np
 import warnings
-from itertools import combinations
+from itertools import combinations, permutations, product
 import scipy.stats
 from scipy.optimize import shgo
 from . import distributions
 from ._continuous_distns import chi2, norm
-from scipy.special import gamma, kv, gammaln
+from scipy.special import gamma, kv, gammaln, comb, factorial
 from . import _wilcoxon_data
+import scipy.stats._bootstrap as _bootstrap
+from scipy._lib._util import check_random_state
 from ._hypotests_pythran import _Q, _P, _a_ij_Aij_Dij2
 
 __all__ = ['epps_singleton_2samp', 'cramervonmises', 'somersd',
-           'barnard_exact', 'boschloo_exact', 'cramervonmises_2samp']
+           'barnard_exact', 'boschloo_exact', 'cramervonmises_2samp',
+           'permutation_test']
 
 Epps_Singleton_2sampResult = namedtuple('Epps_Singleton_2sampResult',
                                         ('statistic', 'pvalue'))
@@ -638,7 +641,7 @@ def _all_partitions(nx, ny):
         y = z[mask]
         yield x, y
 
-        
+
 def _compute_log_combinations(n):
     """Compute all log combination of C(n, k)."""
     gammaln_arr = gammaln(np.arange(n + 1) + 1)
@@ -984,8 +987,8 @@ def boschloo_exact(table, alternative="two-sided", n=32):
     - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 \neq p_2`,
       with `alternative` = "two-sided" (default one)
 
-    Boschloo's exact test uses the p-value of Fisher's exact test as a 
-    statistic, and Boschloo's p-value is the probability under the null 
+    Boschloo's exact test uses the p-value of Fisher's exact test as a
+    statistic, and Boschloo's p-value is the probability under the null
     hypothesis of observing such an extreme value of this statistic.
 
     Boschloo's and Barnard's are both more powerful than Fisher's exact
@@ -1393,3 +1396,692 @@ def cramervonmises_2samp(x, y, method='auto'):
             p = max(0, 1. - _cdf_cvm_inf(tn))
 
     return CramerVonMisesResult(statistic=t, pvalue=p)
+
+
+attributes = ('statistic', 'pvalue', 'null_distribution')
+PermutationTestResult = make_dataclass('PermutationTestResult', attributes)
+
+
+# Functions like these are under review in gh-13312; nice to have them here
+# Remove when they are added to scipy/stats/stats.py or a utility file
+def _broadcast_arrays(arrays, axis=None):
+    """
+    Broadcast shapes of arrays, ignoring incompatibility of specified axes
+    """
+    new_shapes = _broadcast_array_shapes(arrays, axis=axis)
+    if axis is None:
+        new_shapes = [new_shapes]*len(arrays)
+    return [np.broadcast_to(array, new_shape)
+            for array, new_shape in zip(arrays, new_shapes)]
+
+
+def _broadcast_array_shapes(arrays, axis=None):
+    """
+    Broadcast shapes of arrays, ignoring incompatibility of specified axes
+    """
+    shapes = [np.asarray(arr).shape for arr in arrays]
+    return _broadcast_shapes(shapes, axis)
+
+
+def _broadcast_shapes(shapes, axis=None):
+    """
+    Broadcast shapes, ignoring incompatibility of specified axes
+    """
+    # First, ensure all shapes have same number of dimensions by prepending 1s.
+    n_dims = max([len(shape) for shape in shapes])
+    new_shapes = np.ones((len(shapes), n_dims), dtype=int)
+    for row, shape in zip(new_shapes, shapes):
+        row[len(row)-len(shape):] = shape  # can't use negative indices (-0:)
+
+    # Remove the shape elements of the axes to be ignored, but remember them.
+    if axis is not None:
+        axis = np.atleast_1d(axis)
+        axis[axis < 0] = n_dims + axis[axis < 0]
+        removed_shapes = new_shapes[:, axis]
+        new_shapes = np.delete(new_shapes, axis, axis=1)
+
+    # If arrays are broadcastable, shape elements that are 1 may be replaced
+    # with a corresponding non-1 shape element. Assuming arrays are
+    # broadcastable, that final shape element can be found with:
+    new_shape = np.max(new_shapes, axis=0)
+    # except in case of an empty array:
+    new_shape *= new_shapes.all(axis=0)
+
+    # Among all arrays, there can only be one unique non-1 shape element.
+    # Therefore, if any non-1 shape element does not match what we found
+    # above, the arrays must not be broadcastable after all.
+    if np.any(~((new_shapes == 1) | (new_shapes == new_shape))):
+        raise ValueError("Array shapes are incompatible for broadcasting.")
+
+    if axis is not None:
+        # Add back the shape elements that were ignored
+        new_axis = axis - np.arange(len(axis))
+        new_shapes = [tuple(np.insert(new_shape, new_axis, removed_shape))
+                      for removed_shape in removed_shapes]
+        return new_shapes
+    else:
+        return tuple(new_shape)
+
+
+def _all_partitions_concatenated(ns):
+    """
+    Generate all partitions of indices of groups of given sizes, concatenated
+    """
+    def all_partitions(z, n):
+        for c in combinations(z, n):
+            x0 = set(c)
+            x1 = z - x0
+            yield [x0, x1]
+
+    def all_partitions_n(z, ns):
+        if len(ns) == 0:
+            yield [z]
+            return
+        for c in all_partitions(z, ns[0]):
+            for d in all_partitions_n(c[1], ns[1:]):
+                yield c[0:1] + d
+
+    z = set(range(np.sum(ns)))
+    for partitioning in all_partitions_n(z, ns[:]):
+        x = np.concatenate([list(partition)
+                            for partition in partitioning]).astype(int)
+        yield x
+
+
+def _batch_generator(iterable, batch):
+    """A generator that yields batches of elements from an iterable"""
+    iterator = iter(iterable)
+    if batch <= 0:
+        raise ValueError("`batch` must be positive.")
+    while True:
+        z = []
+        try:
+            # get elements from iterator `batch` at a time
+            for i in range(batch):
+                z.append(next(iterator))
+            yield z
+        except StopIteration:
+            # when there are no more elements, yield the final batch and stop
+            if z:
+                yield z
+            break
+
+
+def _calculate_null_both(data, statistic, n_permutations, batch,
+                         random_state=None):
+    """
+    Calculate null distribution for independent sample tests.
+    """
+    n_samples = len(data)
+
+    # compute number of permutations
+    # (distinct partitions of data into samples of these sizes)
+    n_obs_i = [sample.shape[-1] for sample in data]  # observations per sample
+    n_obs_ic = np.cumsum(n_obs_i)
+    n_obs = n_obs_ic[-1]  # total number of observations
+    n_max = np.prod([comb(n_obs_ic[i], n_obs_ic[i-1])
+                     for i in range(n_samples-1, 0, -1)])
+
+    _memory_check(data, min(n_permutations, n_max))
+
+    # perm_generator is an iterator that produces permutations of indices
+    # from 0 to n_obs. We'll concatenate the samples, use these indices to
+    # permute the data, then split the samples apart again.
+    if n_permutations >= n_max:
+        n_permutations = n_max
+        perm_generator = _all_partitions_concatenated(n_obs_i)
+    else:
+        # Neither RandomState.permutation nor Generator.permutation
+        # can permute axis-slices independently. If this feature is
+        # added in the future, batches of the desired size should be
+        # generated in a single call.
+        perm_generator = (random_state.permutation(n_obs)
+                          for i in range(n_permutations))
+
+    batch = batch or int(n_permutations)
+    null_distribution = []
+
+    # First, concatenate all the samples. In batches, permute samples with
+    # indices produced by the `perm_generator`, split them into new samples of
+    # the original sizes, compute the statistic for each batch, and add these
+    # statistic values to the null distribution.
+    data = np.concatenate(data, axis=-1)
+    for indices in _batch_generator(perm_generator, batch=batch):
+        indices = np.array(indices)
+
+        # `indices` is 2D: each row is a permutation of the indices.
+        # We use it to index `data` along its last axis, which corresponds
+        # with observations.
+        # After indexing, the second to last axis of `data_batch` corresponds
+        # with permutations, and the last axis corresponds with observations.
+        data_batch = data[..., indices]
+
+        # Move the permutation axis to the front: we'll concatenate a list
+        # of batched statistic values along this zeroth axis to form the
+        # null distribution.
+        data_batch = np.moveaxis(data_batch, -2, 0)
+        data_batch = np.split(data_batch, n_obs_ic[:-1], axis=-1)
+        null_distribution.append(statistic(*data_batch, axis=-1))
+    null_distribution = np.concatenate(null_distribution, axis=0)
+
+    return null_distribution, n_permutations
+
+
+def _calculate_null_pairings(data, statistic, n_permutations, batch,
+                             random_state=None):
+    """
+    Calculate null distribution for association tests.
+    """
+    n_samples = len(data)
+
+    # compute number of permutations (factorial(n) permutations of each sample)
+    n_obs_sample = data[0].shape[-1]  # observations per sample; same for each
+    n_max = factorial(n_obs_sample)**n_samples
+
+    _memory_check(data, min(n_permutations, n_max))
+
+    # `perm_generator` is an iterator that produces a list of permutations of
+    # indices from 0 to n_obs_sample, one for each sample.
+    if n_permutations >= n_max:
+        n_permutations = n_max
+        # cartesian product of the sets of all permutations of indices
+        perm_generator = product(*(permutations(range(n_obs_sample))
+                                   for i in range(n_samples)))
+    else:
+        # Separate random permutations of indices for each sample.
+        # Again, it would be nice if RandomState/Generator.permutation
+        # could permute each axis-slice separately.
+        perm_generator = ([random_state.permutation(n_obs_sample)
+                           for i in range(n_samples)]
+                          for j in range(n_permutations))
+
+    batch = batch or int(n_permutations)
+    null_distribution = []
+
+    for indices in _batch_generator(perm_generator, batch=batch):
+        indices = np.array(indices)
+
+        # `indices` is 3D: the zeroth axis is for permutations, the next is
+        # for samples, and the last is for observations. Swap the first two
+        # to make the zeroth axis correspond with samples, as it does for
+        # `data`.
+        indices = np.swapaxes(indices, 0, 1)
+
+        # When we're done, `data_batch` will be a list of length `n_samples`.
+        # Each element will be a batch of random permutations of one sample.
+        # The zeroth axis of each batch will correspond with permutations,
+        # and the last will correspond with observations. (This makes it
+        # easy to pass into `statistic`.)
+        data_batch = [None]*n_samples
+        for i in range(n_samples):
+            data_batch[i] = data[i][..., indices[i]]
+            data_batch[i] = np.moveaxis(data_batch[i], -2, 0)
+
+        null_distribution.append(statistic(*data_batch, axis=-1))
+    null_distribution = np.concatenate(null_distribution, axis=0)
+
+    return null_distribution, n_permutations
+
+
+def _calculate_null_samples(data, statistic, n_permutations, batch,
+                            random_state=None):
+    """
+    Calculate null distribution for paired-sample tests.
+    """
+    n_samples = len(data)
+
+    # By convention, the meaning of the "samples" permutations type for
+    # data with only one sample is to flip the sign of the observations.
+    # Achieve this by adding a second sample - the negative of the original.
+    if n_samples == 1:
+        data = [data[0], -data[0]]
+
+    # The "samples" permutation strategy is the same as the "pairings"
+    # strategy except the roles of samples and observations are flipped.
+    # So swap these axes, then we'll use the function for the "pairings"
+    # strategy to do all the work!
+    data = np.swapaxes(data, 0, -1)
+
+    # (Of course, the user's statistic doesn't know what we've done here,
+    # so we need to pass it what it's expecting.)
+    def statistic_wrapped(*data, axis):
+        data = np.swapaxes(data, 0, -1)
+        if n_samples == 1:
+            data = data[0:1]
+        return statistic(*data, axis=axis)
+
+    return _calculate_null_pairings(data, statistic_wrapped, n_permutations,
+                                    batch, random_state)
+
+
+def _permutation_test_iv(data, statistic, permutation_type, vectorized,
+                         permutations, batch, alternative, axis, random_state):
+    """Input validation for `permutation_test`."""
+
+    axis_int = int(axis)
+    if axis != axis_int:
+        raise ValueError("`axis` must be an integer.")
+
+    permutation_types = {'samples', 'pairings', 'both'}
+    permutation_type = permutation_type.lower()
+    if permutation_type not in permutation_types:
+        raise ValueError(f"`permutation_type` must be in {permutation_types}.")
+
+    if vectorized not in {True, False}:
+        raise ValueError("`vectorized` must be `True` or `False`.")
+
+    message = "`data` must be a tuple containing at least two samples"
+    try:
+        if len(data) < 2 and permutation_type == "both":
+            raise ValueError(message)
+    except TypeError:
+        raise TypeError(message)
+
+    data = _broadcast_arrays(data, axis)
+    data_iv = []
+    for sample in data:
+        sample = np.atleast_1d(sample)
+        if sample.shape[axis] <= 1:
+            raise ValueError("each sample in `data` must contain two or more "
+                             "observations along `axis`.")
+        sample = np.moveaxis(sample, axis_int, -1)
+        data_iv.append(sample)
+
+    permutations_int = (int(permutations) if not np.isinf(permutations)
+                        else np.inf)
+    if permutations != permutations_int or permutations_int <= 0:
+        raise ValueError("`permutations` must be a positive integer.")
+
+    if batch is None:
+        batch_iv = batch
+    else:
+        batch_iv = int(batch)
+        if batch != batch_iv or batch_iv <= 0:
+            raise ValueError("`batch` must be a positive integer or None.")
+
+    alternatives = {'two-sided', 'greater', 'less'}
+    alternative = alternative.lower()
+    if alternative not in alternatives:
+        raise ValueError(f"`alternative` must be in {alternatives}")
+
+    random_state = check_random_state(random_state)
+
+    return (data_iv, statistic, permutation_type, vectorized, permutations_int,
+            batch_iv, alternative, axis_int, random_state)
+
+
+# This is mostly to protect me while developing. Eventually we want to batch
+# the computations, so this will be less of an issue, but if it sounds like a
+# good feature, let me know.
+def _memory_check(data, permutations):
+    needed = sum([sample.size for sample in data]) * permutations * 8
+    available = 1e9  # psutil.virtual_memory()[1] # can't use in SciPy
+    if needed > 0.8 * available:
+        raise MemoryError("Not enough memory available. "
+                          "Consider reducing the batch size.")
+
+
+def permutation_test(data, statistic, *, permutation_type='both',
+                     vectorized=False, permutations=np.inf, batch=None,
+                     alternative="two-sided", axis=0, random_state=None):
+    r"""
+    Performs a permutation test of a given statistic on provided data.
+
+    For independent sample statistics, the null hypothesis is that the data are
+    randomly sampled from the same distribution.
+    For paired sample statistics, two null hypothesis can be tested:
+    that the data are paired at random or that the data are assigned to samples
+    at random.
+
+    Parameters
+    ----------
+    data : iterable of array-like
+        Contains the samples, each of which is an array of observations.
+        Dimensions of sample arrays must be compatible for broadcasting except
+        along `axis`.
+    statistic : callable
+        Statistic for which the p-value of the hypothesis test is to be
+        calculated. `statistic` must be a callable that accepts samples
+        as separate arguments (e.g. ``statistic(*data)``) and returns the
+        resulting statistic.
+        If `vectorized` is set ``True``, `statistic` must also accept a keyword
+        argument `axis` and be vectorized to compute the statistic along the
+        provided `axis` of the sample arrays.
+    permutation_type : {'both', 'samples', 'pairings'}, optional
+        The type of permutations to be performed, in accordance with the
+        null hypothesis. The first two permutation types are for paired sample
+        statistics, in which all samples contain the same number of
+        observations and observations with corresponding indices along `axis`
+        are considered to be paired; the third is for independent sample
+        statistics.
+
+        - ``'samples'`` : observations are assigned to different samples
+          but remain paired with the same observations from other samples.
+          This permutation type is appropriate for paired sample hypothesis
+          tests such as the Wilcoxon signed-rank test and the paired t-test.
+        - ``'pairings'`` : observations are paired with different observations,
+          but they remain within the same sample. This permutation type is
+          appropriate for association/correlation tests with statistics such
+          as Spearman's :math:`\rho`, Kendall's :math:`\tau`, and Pearson's
+          :math:`r`.
+        - ``'both'`` (default) : observations are assigned to different samples
+          without preserving pairs. Samples may contain different numbers of
+          observations. This permutation type is appropriate for
+          independent sample hypothesis tests such as the Mann-Whitney
+          :math:`U` test and the independent sample t-test.
+
+          Please see the Notes section below for more detailed descriptions
+          of the permutation types.
+
+    vectorized : bool, default: ``False``
+        By default, `statistic` is assumed to calculate the statistic only for
+        1D arrays contained in `data`. If `vectorized` is set ``True``,
+        `statistic` must also accept a keyword argument `axis` and be
+        vectorized to compute the statistic along the provided `axis` of the ND
+        arrays in `data`. Use of a vectorized statistic can reduce computation
+        time.
+    permutations : int, default: ``np.inf``
+        Number of random permutations used to approximate the null
+        distribution. If greater than or equal to the number of distinct
+        permutations, the exact null distribution will be computed.
+        Note that the number of distinct permutations grows very rapidly with
+        the sizes of samples, so exact tests are feasible only for very small
+        data sets.
+    batch : int, optional
+        The number of permutations to process in each call to `statistic`.
+        Memory usage is O(`batch`*``n``), where ``n`` is the total size
+        of all samples, regardless of the value of `vectorized`. Default is
+        ``None``, in which case ``batch`` is the number of permutations.
+    alternative : {'two-sided', 'less', 'greater'}, optional
+        The alternative hypothesis for which the p-value is calculated.
+        For each alternative, the p-value is defined as follows.
+
+        - ``'greater'`` : the percentage of the null distribution that is
+          greater than or equal to the observed value of the test statistic.
+        - ``'less'`` : the percentage of the null distribution that is
+          less than or equal to the observed value of the test statistic.
+        - ``'two-sided'`` (default) : twice the smaller of the p-values above.
+
+    axis : int, default: 0
+        The axis of the (broadcasted) samples over which to calculate the
+        statistic. If samples have a different number of dimensions,
+        singleton dimensions are prepended to samples with fewer dimensions
+        before `axis` is considered.
+    random_state : {None, int, `numpy.random.Generator`,
+                    `numpy.random.RandomState`}, optional
+
+        Pseudorandom number generator state used to generate resamples.
+
+        If `random_state` is ``None`` (default), the
+        `numpy.random.RandomState` singleton is used.
+        If `random_state` is an int, a new ``RandomState`` instance is used,
+        seeded with `random_state`.
+        If `random_state` is already a ``Generator`` or ``RandomState``
+        instance then that instance is used.
+
+    Returns
+    -------
+    statistic : float or ndarray
+        The observed test statistic of the data.
+    pvalue : float or ndarray
+        The p-value for the given alternative.
+    null_distribution : ndarray
+        The values of the test statistic generated under the null hypothesis.
+
+    Notes
+    -----
+
+    The three types of permutation tests supported by this function are
+    described below.
+
+    **Unpaired statistics** (``permutation_type='both'``):
+
+    The null hypothesis associated with this permutation type is that all
+    observations are sampled from the same underlying distribution and that
+    they have been assigned to one of the samples at random.
+
+    Suppose ``data`` contains two samples; e.g. ``a, b = data``.
+    When ``1 < permutations < binom(n, k)``, where
+
+    * ``k`` is the number of observations in ``a``,
+    * ``n`` is the total number of observations in ``a`` and ``b``, and
+    * ``binom(n, k)`` is the binomial coefficient (``n`` choose ``k``),
+
+    the data are pooled (concatenated), randomly assigned to either the first
+    or second sample, and the statistic is calculated. This process is
+    performed repeatedly, `permutation` times, generating a distribution of the
+    statistic under the null hypothesis. The statistic of the original
+    data is compared to this distribution to determine the p-value.
+
+    When ``permutations >= binom(n, k)``, an exact test is performed: the data
+    are *partitioned* between the samples in each distinct way exactly once,
+    and the exact null distribution is formed.
+    Note that for a given partitioning of the data between the samples,
+    only one ordering/permutation of the data *within* each sample is
+    considered. For statistics that do not depend on the order of the data
+    within samples, this dramatically reduces computational cost without
+    affecting the shape of the null distribution (because the frequency/count
+    of each value is affected by the same factor).
+
+    For ``a = [a1, a2, a3, a4]`` and ``b = [b1, b2, b3]``, an example of this
+    permutation type is ``x = [b3, a1, a2, b2]`` and ``y = [a4, b1, a3]``.
+
+    ``permutation_type='both'`` does not support one-sample statistics, but it
+    can be applied to statistics with more than two samples. In this case, if
+    ``n`` is an array of the number of observations within each sample, the
+    number of distinct partitions is::
+
+        np.product([binom(sum(n[i:]), sum(n[i+1:])) for i in range(len(n)-1)])
+
+    **Paired statistics, permute pairings** (``permutation_type='pairings'``):
+
+    The null hypothesis associated with this permutation type is that
+    observations within each sample are drawn from the same underlying
+    distribution and that pairings with elements of other samples are
+    assigned at random.
+
+    Suppose ``data`` contains only one sample; e.g. ``a, = data``, and we
+    wish to consider all possible pairings of elements of ``a`` with elements
+    of a second sample, ``b``. Let ``n`` be the number of observations in
+    ``a``, which must also equal the number of observations in ``b``.
+
+    When ``1 < permutations < factorial(n)``, the elements of ``a`` are
+    randomly permuted. The user-supplied statistic accepts one data argument,
+    say ``a_perm``, and calculates the statistic considering ``a_perm`` and
+    ``b``. This process is performed repeatedly, `permutation` times,
+    generating a distribution of the statistic under the null hypothesis.
+    The statistic of the original data is compared to this distribution to
+    determine the p-value.
+
+    When ``permutations >= factorial(n)``, an exact test is performed:
+    ``a`` is permuted in each distinct way exactly once. Therefore, the
+    `statistic` is computed for each unique pairing of samples between ``a``
+    and ``b`` exactly once.
+
+    For ``a = [a1, a2, a3]`` and ``b = [b1, b2, b3]``, an example of this
+    permutation type is ``a_perm = [a3, a1, a2]`` while ``b`` is left
+    in its original order.
+
+    ``permutation_type='pairings'`` supports ``data`` containing any number
+    of samples, each of which must contain the same number of observations.
+    All samples provided in ``data`` are permuted *independently*. Therefore,
+    if ``m`` is the number of samples and ``n`` is the number of observations
+    within each sample, then the number of permutations in an exact test is::
+
+        factorial(n)**m
+
+    Note that if a two-sample statistic, for example, does not inherently
+    depend on the order in which observations are provided - only on the
+    *pairings* of observations - then only one of the two samples should be
+    provided in ``data``. This dramatically reduces computational cost without
+    affecting the shape of the null distribution (because the frequency/count
+    of each value is affected by the same factor).
+
+    **Paired statistics, permute samples** (``permutation_type='samples'``):
+
+    The null hypothesis associated with this permutation type is that
+    observations within each pair are drawn from the same underlying
+    distribution and that the sample to which they are assigned is random.
+
+    Suppose ``data`` contains two samples; e.g. ``a, b = data``.
+    Let ``n`` be the number of observations in ``a``, which must also equal
+    the number of observations in ``b``.
+
+    When ``1 < permutations < 2**n``, the elements of ``a`` are ``b`` are
+    randomly swapped between samples (maintaining their pairings) and the
+    statistic is calculated. This process is performed repeatedly,
+    `permutation` times,  generating a distribution of the statistic under the
+    null hypothesis. The statistic of the original data is compared to this
+    distribution to determine the p-value.
+
+    When ``permutations >= 2**n``, an exact test is performed: the observations
+    are assigned to the two samples in each distinct way (while maintaining
+    pairings) exactly once.
+
+    For ``a = [a1, a2, a3]`` and ``b = [b1, b2, b3]``, an example of this
+    permutation type is ``x = [b1, a2, b3]`` and ``y = [a1, b2, a3]``.
+
+    ``permutation_type='samples'`` supports ``data`` containing any number
+    of samples, each of which must contain the same number of observations.
+    If ``data`` contains more than one sample, paired observations within
+    ``data`` are exchanged between samples *independently*. Therefore, if ``m``
+    is the number of samples and ``n`` is the number of observations within
+    each sample, then the number of permutations in an exact test is::
+
+        factorial(m)**n
+
+    Several paired-sample statistical tests, such as the Wilcoxon signed rank
+    test and paired-sample t-test, can be performed considering only the
+    *difference* between two paired elements. Accordingly, if ``data`` contains
+    only one sample, then the null distribution is formed by independently
+    changing the *sign* of each observation.
+
+    Examples
+    --------
+
+    Suppose we wish to test whether two samples are drawn from the same
+    distribution. Assume that the underlying distributions are unknown to us,
+    and that before observing the data, we hypothesized that the mean of the
+    first sample would be less than that of the second sample. We decide that
+    we will use the difference between the sample means as a test statistic,
+    and we will consider a p-value of 0.05 to be statistically significant.
+
+    For efficiency, we write the function defining the test statistic in a
+    vectorized fashion: the samples ``x`` and ``y`` can be ND arrays, and the
+    statistic will be calculated for each axis-slice along `axis`.
+
+    >>> def statistic(x, y, axis):
+    ...     return np.mean(x, axis=axis) - np.mean(y, axis=axis)
+
+    After collecting our data, we calculate the observed value of the test
+    statistic.
+
+    >>> from scipy.stats import norm
+    >>> rng = np.random.default_rng()
+    >>> x = norm.rvs(size=5, random_state=rng)
+    >>> y = norm.rvs(size=6, loc = 3, random_state=rng)
+    >>> statistic(x, y, 0)
+    -3.5411688580987266
+
+    Indeed, the test statistic is negative, suggesting that the true mean of
+    the distribution underlying ``x`` is less than that of the distribution
+    underlying ``y``. To determine the probability of this occuring by chance
+    if the two samples were drawn from the same distribution, we perform
+    a permutation test.
+
+    >>> from scipy.stats import permutation_test
+    >>> # because our statistic is vectorized, we pass `vectorized=True`
+    >>> res = permutation_test((x, y), statistic, vectorized=True,
+    ...                        alternative='less')
+    >>> print(res.statistic)
+    -3.5411688580987266
+    >>> print(res.pvalue)
+    0.004329004329004329
+
+    The probability of obtaining a test statistic less than or equal to the
+    observed value under the null hypothesis is 0.4329%. This is less than our
+    chosen threshold of 5%, so we consider this to to be significant evidence
+    against the null hypothesis in favor of the alternative.
+
+    Because the size of the samples above was small, `permutation_test` could
+    perform an exact test. For larger samples, we resort to a randomized
+    permutation test.
+
+    >>> x = norm.rvs(size=100, random_state=rng)
+    >>> y = norm.rvs(size=120, loc=0.3, random_state=rng)
+    >>> res = permutation_test((x, y), statistic, permutations=100000,
+    ...                        vectorized=True, alternative='less',
+    ...                        random_state=rng)
+    >>> print(res.statistic)
+    -0.5230459671240913
+    >>> print(res.pvalue)
+    0.00016
+
+    The approximate probability of obtaining a test statistic less than or
+    equal to the observed value under the null hypothesis is 0.0225%. This is
+    again less than our chosen threshold of 5%, so again we have significant
+    evidence to reject the null hypothesis in favor of the alternative.
+
+    For large samples and number of permutations, the result is comparable to
+    that of the corresponding asymptotic test, the independent sample t-test.
+
+    >>> from scipy.stats import ttest_ind
+    >>> res_asymptotic = ttest_ind(x, y, alternative='less')
+    >>> print(res_asymptotic.pvalue)
+    0.00012688101537979522
+
+    The permutation distribution of the test statistic is provided for
+    further investigation.
+
+    >>> import matplotlib.pyplot as plt
+    >>> plt.hist(res.null_distribution, bins=50)
+    >>> plt.title("Permutation distribution of test statistic")
+    >>> plt.xlabel("Value of Statistic")
+    >>> plt.ylabel("Frequency")
+
+    """
+    args = _permutation_test_iv(data, statistic, permutation_type, vectorized,
+                                permutations, batch, alternative, axis,
+                                random_state)
+    (data, statistic, permutation_type, vectorized, permutations, batch,
+     alternative, axis, random_state) = args
+
+    if not vectorized:
+        statistic_vectorized = _bootstrap._vectorize_statistic(statistic)
+    else:
+        statistic_vectorized = statistic
+
+    observed = statistic_vectorized(*data, axis=-1)
+
+    null_calculators = {"pairings": _calculate_null_pairings,
+                        "samples": _calculate_null_samples,
+                        "both": _calculate_null_both}
+    null_calculator_args = (data, statistic_vectorized, permutations,
+                            batch, random_state)
+    calculate_null = null_calculators[permutation_type]
+    null_distribution, permutations = calculate_null(*null_calculator_args)
+
+    def less(null_distribution, observed):
+        cmps = null_distribution <= observed
+        pvalues = cmps.sum(axis=0) / permutations
+        return pvalues
+
+    def greater(null_distribution, observed):
+        cmps = null_distribution >= observed
+        pvalues = cmps.sum(axis=0) / permutations
+        return pvalues
+
+    def two_sided(null_distribution, observed):
+        pvalues_less = less(null_distribution, observed)
+        pvalues_greater = greater(null_distribution, observed)
+        pvalues = np.minimum(pvalues_less, pvalues_greater) * 2
+        return pvalues
+
+    compare = {"less": less,
+               "greater": greater,
+               "two-sided": two_sided}
+
+    pvalues = compare[alternative](null_distribution, observed)
+    pvalues = np.clip(pvalues, 0, 1)
+
+    return PermutationTestResult(observed, pvalues, null_distribution)
diff --git a/scipy/stats/tests/test_hypotests.py b/scipy/stats/tests/test_hypotests.py
index d4732a3ffc6d..234e115c3d11 100644
--- a/scipy/stats/tests/test_hypotests.py
+++ b/scipy/stats/tests/test_hypotests.py
@@ -13,9 +13,12 @@
 from scipy.stats._hypotests import (epps_singleton_2samp, cramervonmises,
                                     _cdf_cvm, cramervonmises_2samp,
                                     _pval_cvm_2samp_exact, barnard_exact,
-                                    boschloo_exact)
+                                    boschloo_exact, permutation_test,
+                                    _all_partitions_concatenated)
 from scipy.stats._mannwhitneyu import mannwhitneyu, _mwu_state
+import scipy.stats._bootstrap as _bootstrap
 from .common_tests import check_named_results
+from scipy import special
 
 
 class TestEppsSingleton:
@@ -1010,6 +1013,7 @@ class TestBoschlooExact:
     """Some tests to show that boschloo_exact() works correctly."""
 
     ATOL = 1e-7
+
     @pytest.mark.parametrize(
         "input_sample,expected",
         [
@@ -1138,6 +1142,7 @@ def test_row_or_col_zero(self, input_sample, expected):
         assert_equal(pvalue, expected[0])
         assert_equal(statistic, expected[1])
 
+
 class TestCvm_2samp:
     def test_invalid_input(self):
         x = np.arange(10).reshape((2, 5))
@@ -1226,3 +1231,629 @@ def test_same_input(self):
         # check exact p-value
         res = cramervonmises_2samp(x[:4], x[:4])
         assert_equal((res.statistic, res.pvalue), (0.0, 1.0))
+
+
+def test_all_partitions_concatenated():
+    # make sure that _all_paritions_concatenated produces the correct number
+    # of partitions of the data into samples of the given sizes and that
+    # all are unique
+    n = np.array([3, 2, 4], dtype=int)
+    nc = np.cumsum(n)
+
+    all_partitions = set()
+    counter = 0
+    for partition_concatenated in _all_partitions_concatenated(n):
+        counter += 1
+        partitioning = np.split(partition_concatenated, nc[:-1])
+        all_partitions.add(tuple([frozenset(i) for i in partitioning]))
+
+    expected = np.product([special.binom(sum(n[i:]), sum(n[i+1:]))
+                           for i in range(len(n)-1)])
+
+    assert_equal(counter, expected)
+    assert_equal(len(all_partitions), expected)
+
+
+class TestPermutationTest:
+
+    rtol = 1e-14
+
+    # -- Input validation -- #
+
+    def test_permutation_test_iv(self):
+
+        def stat(x, y, axis):
+            return stats.ttest_ind((x, y), axis).statistic
+
+        message = "each sample in `data` must contain two or more ..."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1]), stat)
+
+        message = "`data` must be a tuple containing at least two samples"
+        with pytest.raises(ValueError, match=message):
+            permutation_test((1,), stat)
+        with pytest.raises(TypeError, match=message):
+            permutation_test(1, stat)
+
+        message = "`axis` must be an integer."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, axis=1.5)
+
+        message = "`permutation_type` must be in..."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat,
+                             permutation_type="ekki")
+
+        message = "`vectorized` must be `True` or `False`."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, vectorized=1.5)
+
+        message = "`permutations` must be a positive integer."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, permutations=-1000)
+
+        message = "`permutations` must be a positive integer."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, permutations=1000.5)
+
+        message = "`batch` must be a positive integer or None."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, batch=-1000)
+
+        message = "`batch` must be a positive integer or None."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, batch=1000.5)
+
+        message = "`alternative` must be in..."
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat, alternative='ekki')
+
+        message = "'herring' cannot be used to seed a"
+        with pytest.raises(ValueError, match=message):
+            permutation_test(([1, 2, 3], [1, 2, 3]), stat,
+                             random_state='herring')
+
+    # -- Test Parameters -- #
+    @pytest.mark.parametrize('permutation_type',
+                             ['pairings', 'samples', 'both'])
+    def test_batch(self, permutation_type):
+        # make sure that the `batch` parameter is respected by checking the
+        # maximum batch size provided in calls to `statistic`
+        np.random.seed(0)
+        x = np.random.rand(10)
+        y = np.random.rand(10)
+
+        def statistic(x, y, axis):
+            batch_size = 1 if x.ndim == 1 else len(x)
+            statistic.batch_size = max(batch_size, statistic.batch_size)
+            statistic.counter += 1
+            return np.mean(x, axis=axis) - np.mean(y, axis=axis)
+        statistic.counter = 0
+        statistic.batch_size = 0
+
+        kwds = {'permutations': 1000, 'permutation_type': permutation_type,
+                'vectorized': True, 'random_state': 0}
+        res1 = stats.permutation_test((x, y), statistic, batch=1, **kwds)
+        assert_equal(statistic.counter, 1001)
+        assert_equal(statistic.batch_size, 1)
+
+        statistic.counter = 0
+        res2 = stats.permutation_test((x, y), statistic, batch=50, **kwds)
+        assert_equal(statistic.counter, 21)
+        assert_equal(statistic.batch_size, 50)
+
+        statistic.counter = 0
+        res3 = stats.permutation_test((x, y), statistic, batch=1000, **kwds)
+        assert_equal(statistic.counter, 2)
+        assert_equal(statistic.batch_size, 1000)
+
+        assert_equal(res1.pvalue, res3.pvalue)
+        assert_equal(res2.pvalue, res3.pvalue)
+
+    @pytest.mark.parametrize('permutation_type, exact_size',
+                             [('pairings', special.factorial(3)**2),
+                              ('samples', 2**3),
+                              ('both', special.binom(6, 3))])
+    def test_permutations(self, permutation_type, exact_size):
+        # make sure that the `permutations` parameter is respected by checking
+        # the size of the null distribution
+        np.random.seed(0)
+        x = np.random.rand(3)
+        y = np.random.rand(3)
+
+        def statistic(x, y, axis):
+            return np.mean(x, axis=axis) - np.mean(y, axis=axis)
+
+        kwds = {'permutation_type': permutation_type,
+                'vectorized': True, 'random_state': 0}
+        res = stats.permutation_test((x, y), statistic, permutations=3, **kwds)
+        assert_equal(res.null_distribution.size, 3)
+
+        res = stats.permutation_test((x, y), statistic, **kwds)
+        assert_equal(res.null_distribution.size, exact_size)
+
+    # -- Randomized Permutation Tests -- #
+
+    # To get reasonable accuracy, these next three tests are somewhat slow.
+    # Originally, I had them passing for all combinations of permutation type,
+    # alternative, and RNG, but that takes too long for CI. Instead, split
+    # into three tests, each testing a particular combination of the three
+    # parameters.
+
+    def test_randomized_test_against_exact_both(self):
+        # check that the randomized and exact tests agree to reasonable
+        # precision for permutation_type='both
+
+        alternative, rng = 'less', 0
+
+        nx, ny, permutations = 8, 9, 24000
+        assert special.binom(nx + ny, nx) > permutations
+
+        x = stats.norm.rvs(size=nx)
+        y = stats.norm.rvs(size=ny)
+        data = x, y
+
+        def statistic(x, y, axis):
+            return np.mean(x, axis=axis) - np.mean(y, axis=axis)
+
+        kwds = {'vectorized': True, 'permutation_type': 'both',
+                'batch': 100, 'alternative': alternative, 'random_state': rng}
+        res = permutation_test(data, statistic, **kwds,
+                               permutations=permutations)
+        res2 = permutation_test(data, statistic, **kwds)
+
+        assert res.statistic == res2.statistic
+        assert_allclose(res.pvalue, res2.pvalue, atol=1e-2)
+
+    def test_randomized_test_against_exact_samples(self):
+        # check that the randomized and exact tests agree to reasonable
+        # precision for permutation_type='samples'
+
+        alternative, rng = 'greater', None
+
+        nx, ny, permutations = 15, 15, 32000
+        assert 2**nx > permutations
+
+        x = stats.norm.rvs(size=nx)
+        y = stats.norm.rvs(size=ny)
+        data = x, y
+
+        def statistic(x, y, axis):
+            return np.mean(x - y, axis=axis)
+
+        kwds = {'vectorized': True, 'permutation_type': 'samples',
+                'batch': 100, 'alternative': alternative, 'random_state': rng}
+        res = permutation_test(data, statistic, **kwds,
+                               permutations=permutations)
+        res2 = permutation_test(data, statistic, **kwds)
+
+        assert res.statistic == res2.statistic
+        assert_allclose(res.pvalue, res2.pvalue, atol=1e-2)
+
+    def test_randomized_test_against_exact_pairings(self):
+        # check that the randomized and exact tests agree to reasonable
+        # precision for permutation_type='pairings'
+
+        alternative = 'two-sided'
+        try:
+            rng = np.random.default_rng(1)
+        except AttributeError:
+            rng = np.random.RandomState(1)
+
+        nx, ny, permutations = 8, 8, 40000
+        assert special.factorial(nx) > permutations
+
+        x = stats.norm.rvs(size=nx)
+        y = stats.norm.rvs(size=ny)
+        data = [x]
+
+        def statistic1d(x):
+            return stats.pearsonr(x, y)[0]
+
+        statistic = _bootstrap._vectorize_statistic(statistic1d)
+
+        kwds = {'vectorized': True, 'permutation_type': 'samples',
+                'batch': 100, 'alternative': alternative, 'random_state': rng}
+        res = permutation_test(data, statistic, **kwds,
+                               permutations=permutations)
+        res2 = permutation_test(data, statistic, **kwds)
+
+        assert res.statistic == res2.statistic
+        assert_allclose(res.pvalue, res2.pvalue, atol=1e-2)
+
+    @pytest.mark.parametrize('alternative', ('less', 'greater'))
+    # Different conventions for two-sided p-value here VS ttest_ind.
+    # Eventually, we can add multiple options for the two-sided alternative
+    # here in permutation_test.
+    @pytest.mark.parametrize('permutations', (30, 1e9))
+    @pytest.mark.parametrize('axis', (0, 1, 2))
+    def test_against_permutation_ttest(self, alternative, permutations, axis):
+        # check that this function and ttest_ind with permutations give
+        # essentially identical results.
+
+        x = np.arange(3*4*5).reshape(3, 4, 5)
+        y = np.moveaxis(np.arange(4)[:, None, None], 0, axis)
+
+        res1 = stats.ttest_ind(x, y, permutations=permutations, axis=axis,
+                               random_state=0, alternative=alternative)
+
+        def statistic(x, y, axis):
+            return stats.ttest_ind(x, y, axis=axis).statistic
+
+        res2 = permutation_test((x, y), statistic, vectorized=True,
+                                permutations=permutations,
+                                alternative=alternative, axis=axis,
+                                random_state=0)
+
+        assert_allclose(res1.statistic, res2.statistic, rtol=self.rtol)
+        assert_allclose(res1.pvalue, res2.pvalue, rtol=self.rtol)
+
+    # -- Independent (Unpaired) Sample Tests -- #
+
+    def test_against_kstest(self):
+        np.random.seed(0)
+        x = stats.norm.rvs(size=4, scale=1)
+        y = stats.norm.rvs(size=5, loc=3, scale=3)
+
+        alternative = 'greater'
+        expected = stats.ks_2samp(x, y, alternative=alternative, mode='exact')
+
+        def statistic1d(x, y):
+            return stats.ks_2samp(x, y, mode='asymp',
+                                  alternative=alternative).statistic
+
+        res = permutation_test((x, y), statistic1d, alternative=alternative)
+
+        assert_allclose(res.statistic, expected.statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected.pvalue, rtol=self.rtol)
+
+    @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided"))
+    def test_against_ansari(self, alternative):
+        np.random.seed(0)
+        x = stats.norm.rvs(size=4, scale=1)
+        y = stats.norm.rvs(size=5, scale=3)
+
+        # ansari has a different convention for 'alternative'
+        alternative_correspondence = {"less": "greater",
+                                      "greater": "less",
+                                      "two-sided": "two-sided"}
+        alternative_scipy = alternative_correspondence[alternative]
+        expected = stats.ansari(x, y, alternative=alternative_scipy)
+
+        def statistic1d(x, y):
+            return stats.ansari(x, y).statistic
+
+        res = permutation_test((x, y), statistic1d, alternative=alternative)
+
+        assert_allclose(res.statistic, expected.statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected.pvalue, rtol=self.rtol)
+
+    @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided"))
+    def test_against_mannwhitneyu(self, alternative):
+        np.random.seed(0)
+        x = stats.uniform.rvs(size=(3, 5, 2), loc=0)
+        y = stats.uniform.rvs(size=(3, 5, 2), loc=0.05)
+
+        expected = stats.mannwhitneyu(x, y, axis=1, alternative=alternative)
+
+        def statistic(x, y, axis):
+            return stats.mannwhitneyu(x, y, axis=axis).statistic
+
+        res = permutation_test((x, y), statistic, vectorized=True,
+                               alternative=alternative, axis=1)
+
+        assert_allclose(res.statistic, expected.statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected.pvalue, rtol=self.rtol)
+
+    def test_against_cvm(self):
+        np.random.seed(0)
+        x = stats.norm.rvs(size=4, scale=1)
+        y = stats.norm.rvs(size=5, loc=3, scale=3)
+
+        expected = stats.cramervonmises_2samp(x, y, method='exact')
+
+        def statistic1d(x, y):
+            return stats.cramervonmises_2samp(x, y,
+                                              method='asymptotic').statistic
+
+        # cramervonmises_2samp has only one alternative, greater
+        res = permutation_test((x, y), statistic1d, alternative='greater')
+
+        assert_allclose(res.statistic, expected.statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected.pvalue, rtol=self.rtol)
+
+    @pytest.mark.slow()
+    @pytest.mark.parametrize('axis', (-1, 2))
+    def test_vectorized_nsamp_ptype_both(self, axis):
+        # Test that permutation_test with permutation_type='both' works
+        # properly for a 3-sample statistic with nd array samples of different
+        # (but compatible) shapes and ndims. Show that exact permutation test
+        # and random permutation tests approximate SciPy's asymptotic pvalues
+        # and that exact and random permutation test results are even closer
+        # to one another (than they are to the asymptotic results).
+        np.random.seed(0)
+
+        # Three samples, different (but compatible) shapes with different ndims
+        x = np.random.rand(3)
+        y = np.random.rand(1, 3, 2)
+        z = np.random.rand(2, 1, 4)
+        data = (x, y, z)
+
+        # Define the statistic (and pvalue for comparison)
+        def statistic1d(*data):
+            return stats.kruskal(*data).statistic
+
+        def pvalue1d(*data):
+            return stats.kruskal(*data).pvalue
+
+        statistic = _bootstrap._vectorize_statistic(statistic1d)
+        pvalue = _bootstrap._vectorize_statistic(pvalue1d)
+
+        # Calculate the expected results
+        x2 = np.broadcast_to(x, (2, 3, 3))  # broadcast manually because
+        y2 = np.broadcast_to(y, (2, 3, 2))  # _vectorize_statistic doesn't
+        z2 = np.broadcast_to(z, (2, 3, 4))
+        expected_statistic = statistic(x2, y2, z2, axis=axis)
+        expected_pvalue = pvalue(x2, y2, z2, axis=axis)
+
+        # Calculate exact and randomized permutation results
+        kwds = {'vectorized': False, 'axis': axis, 'alternative': 'greater',
+                'permutation_type': 'both', 'random_state': 0}
+        res = permutation_test(data, statistic1d, permutations=np.inf, **kwds)
+        res2 = permutation_test(data, statistic1d, permutations=1000, **kwds)
+
+        # Check results
+        assert_allclose(res.statistic, expected_statistic, rtol=self.rtol)
+        assert_allclose(res.statistic, res2.statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected_pvalue, atol=6e-2)
+        assert_allclose(res.pvalue, res2.pvalue, atol=3e-2)
+
+    # -- Paired-Sample Tests -- #
+
+    @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided"))
+    def test_against_wilcoxon(self, alternative):
+        np.random.seed(0)
+        x = stats.uniform.rvs(size=(3, 6, 2), loc=0)
+        y = stats.uniform.rvs(size=(3, 6, 2), loc=0.05)
+
+        # We'll check both 1- and 2-sample versions of the same test;
+        # we expect identical results to wilcoxon in all cases.
+        def statistic_1samp_1d(z):
+            # 'less' ensures we get the same of two statistics every time
+            return stats.wilcoxon(z, alternative='less').statistic
+
+        def statistic_2samp_1d(x, y):
+            return stats.wilcoxon(x, y, alternative='less').statistic
+
+        def test_1d(x, y):
+            return stats.wilcoxon(x, y, alternative=alternative)
+
+        test = _bootstrap._vectorize_statistic(test_1d)
+
+        expected = test(x, y, axis=1)
+        expected_stat = expected[:, 0, :]
+        expected_p = expected[:, 1, :]
+
+        res1 = permutation_test((x-y,), statistic_1samp_1d, vectorized=False,
+                                permutation_type='samples',
+                                alternative=alternative, axis=1)
+
+        res2 = permutation_test((x, y), statistic_2samp_1d, vectorized=False,
+                                permutation_type='samples',
+                                alternative=alternative, axis=1)
+
+        # `wilcoxon` returns a different statistic with 'two-sided'
+        assert_allclose(res1.statistic, res2.statistic, rtol=self.rtol)
+        if alternative != 'two-sided':
+            assert_allclose(res2.statistic, expected_stat, rtol=self.rtol)
+
+        assert_allclose(res2.pvalue, expected_p, rtol=self.rtol)
+        assert_allclose(res1.pvalue, res2.pvalue, rtol=self.rtol)
+
+    @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided"))
+    def test_against_binomtest(self, alternative):
+        np.random.seed(0)
+        x = np.random.randint(0, 2, size=10)
+        x[x == 0] = -1
+        # More naturally, the test would flip elements between 0 and one.
+        # However, permutation_test will flip the _signs_ of the elements.
+        # So we have to work with +1/-1 instead of 1/0.
+
+        def statistic(x, axis=0):
+            return np.sum(x > 0, axis=axis)
+
+        k, n, p = statistic(x), 10, 0.5
+        expected = stats.binomtest(k, n, p, alternative=alternative)
+
+        res = stats.permutation_test((x,), statistic, vectorized=True,
+                                     permutation_type='samples',
+                                     alternative=alternative)
+        assert_allclose(res.pvalue, expected.pvalue, rtol=self.rtol)
+
+    # -- Exact Association Tests -- #
+
+    def test_against_kendalltau(self):
+        np.random.seed(0)
+        x = stats.norm.rvs(size=6)
+        y = x + stats.norm.rvs(size=6)
+
+        expected = stats.kendalltau(x, y, method='exact')
+
+        def statistic1d(x):
+            return stats.kendalltau(x, y, method='asymptotic').correlation
+
+        # kendalltau currently has only one alternative, two-sided
+        res = permutation_test((x,), statistic1d, permutation_type='pairings')
+
+        assert_allclose(res.statistic, expected.correlation, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected.pvalue, rtol=self.rtol)
+
+    @pytest.mark.parametrize('alternative', ('less', 'greater', 'two-sided'))
+    def test_against_fisher_exact(self, alternative):
+
+        def statistic(x,):
+            return np.sum((x == 1) & (y == 1))
+
+        np.random.seed(0)
+        # x and y are binary random variables with some dependence
+        x = (np.random.rand(7) > 0.6).astype(float)
+        y = (np.random.rand(7) + 0.25*x > 0.6).astype(float)
+        tab = stats.contingency.crosstab(x, y)[1]
+
+        res = permutation_test((x,), statistic, permutation_type='pairings',
+                               alternative=alternative)
+        res2 = stats.fisher_exact(tab, alternative=alternative)
+
+        assert_allclose(res.pvalue, res2[1])
+
+    @pytest.mark.slow()
+    @pytest.mark.parametrize('axis', (-2, 1))
+    def test_vectorized_nsamp_ptype_samples(self, axis):
+        # Test that permutation_test with permutation_type='samples' works
+        # properly for a 3-sample statistic with nd array samples of different
+        # (but compatible) shapes and ndims. Show that exact permutation test
+        # reproduces SciPy's exact pvalue and that random permutation test
+        # approximates it.
+
+        np.random.seed(0)
+
+        x = np.random.rand(2, 4, 3)
+        y = np.random.rand(1, 4, 3)
+        z = np.random.rand(2, 4, 1)
+        x = stats.rankdata(x, axis=axis)
+        y = stats.rankdata(y, axis=axis)
+        z = stats.rankdata(z, axis=axis)
+        y = y[0]  # to check broadcast with different ndim
+        data = (x, y, z)
+
+        def statistic1d(*data):
+            return stats.page_trend_test(data, ranked=True,
+                                         method='asymptotic').statistic
+
+        def pvalue1d(*data):
+            return stats.page_trend_test(data, ranked=True,
+                                         method='exact').pvalue
+
+        statistic = _bootstrap._vectorize_statistic(statistic1d)
+        pvalue = _bootstrap._vectorize_statistic(pvalue1d)
+
+        expected_statistic = statistic(*np.broadcast_arrays(*data), axis=axis)
+        expected_pvalue = pvalue(*np.broadcast_arrays(*data), axis=axis)
+
+        kwds = {'vectorized': False, 'axis': axis, 'alternative': 'greater',
+                'permutation_type': 'pairings', 'random_state': 0}
+        res = permutation_test(data, statistic1d, permutations=np.inf, **kwds)
+        res2 = permutation_test(data, statistic1d, permutations=5000, **kwds)
+
+        assert_allclose(res.statistic, expected_statistic, rtol=self.rtol)
+        assert_allclose(res.statistic, res2.statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected_pvalue, rtol=self.rtol)
+        assert_allclose(res.pvalue, res2.pvalue, atol=3e-2)
+
+    # -- Test Against External References -- #
+
+    tie_case_1 = {'x': [1, 2, 3, 4], 'y': [1.5, 2, 2.5],
+                  'expected_less': 0.2000000000,
+                  'expected_2sided': 0.4,  # 2*expected_less
+                  'expected_Pr_gte_S_mean': 0.3428571429,  # see note below
+                  'expected_statistic': 7.5,
+                  'expected_avg': 9.142857, 'expected_std': 1.40698}
+    tie_case_2 = {'x': [111, 107, 100, 99, 102, 106, 109, 108],
+                  'y': [107, 108, 106, 98, 105, 103, 110, 105, 104],
+                  'expected_less': 0.1555738379,
+                  'expected_2sided': 0.3111476758,
+                  'expected_Pr_gte_S_mean': 0.2969971205,  # see note below
+                  'expected_statistic': 32.5,
+                  'expected_avg': 38.117647, 'expected_std': 5.172124}
+
+    @pytest.mark.slow()  # only the second case is slow, really
+    @pytest.mark.parametrize('case', (tie_case_1, tie_case_2))
+    def test_with_ties(self, case):
+        """
+        Results above from SAS PROC NPAR1WAY, e.g.
+
+        DATA myData;
+        INPUT X Y;
+        CARDS;
+        1 1
+        1 2
+        1 3
+        1 4
+        2 1.5
+        2 2
+        2 2.5
+        ods graphics on;
+        proc npar1way AB data=myData;
+            class X;
+            EXACT;
+        run;
+        ods graphics off;
+
+        Note: SAS provides Pr >= |S-Mean|, which is different from our
+        definition of a two-sided p-value.
+
+        """
+
+        x = case['x']
+        y = case['y']
+
+        expected_statistic = case['expected_statistic']
+        expected_less = case['expected_less']
+        expected_2sided = case['expected_2sided']
+        expected_Pr_gte_S_mean = case['expected_Pr_gte_S_mean']
+        expected_avg = case['expected_avg']
+        expected_std = case['expected_std']
+
+        def statistic1d(x, y):
+            return stats.ansari(x, y).statistic
+
+        with np.testing.suppress_warnings() as sup:
+            sup.filter(UserWarning, "Ties preclude use of exact statistic")
+            res = permutation_test((x, y), statistic1d, alternative='less')
+            res2 = permutation_test((x, y), statistic1d,
+                                    alternative='two-sided')
+
+        assert_allclose(res.statistic, expected_statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected_less, atol=1e-10)
+        assert_allclose(res2.pvalue, expected_2sided, atol=1e-10)
+        assert_allclose(res2.null_distribution.mean(), expected_avg, rtol=1e-6)
+        assert_allclose(res2.null_distribution.std(), expected_std, rtol=1e-6)
+
+        # SAS provides Pr >= |S-Mean|; might as well check against that, too
+        S = res.statistic
+        mean = res.null_distribution.mean()
+        n = len(res.null_distribution)
+        Pr_gte_S_mean = np.sum(np.abs(res.null_distribution-mean)
+                               >= np.abs(S-mean))/n
+        assert_allclose(expected_Pr_gte_S_mean, Pr_gte_S_mean)
+
+    @pytest.mark.parametrize('alternative, expected_pvalue',
+                             (('less', 0.9708333333333),
+                              ('greater', 0.05138888888889),
+                              ('two-sided', 0.1027777777778)))
+    def test_against_spearmanr_in_R(self, alternative, expected_pvalue):
+        """
+        Results above from R cor.test, e.g.
+
+        options(digits=16)
+        x <- c(1.76405235, 0.40015721, 0.97873798,
+               2.2408932, 1.86755799, -0.97727788)
+        y <- c(2.71414076, 0.2488, 0.87551913,
+               2.6514917, 2.01160156, 0.47699563)
+        cor.test(x, y, method = "spearm", alternative = "t")
+        """
+        np.random.seed(0)
+        x = stats.norm.rvs(size=6)
+        y = x + stats.norm.rvs(size=6)
+
+        expected_statistic = 0.7714285714285715
+
+        def statistic1d(x):
+            return stats.spearmanr(x, y).correlation
+
+        res = permutation_test((x,), statistic1d,
+                               permutation_type='pairings',
+                               alternative=alternative)
+
+        assert_allclose(res.statistic, expected_statistic, rtol=self.rtol)
+        assert_allclose(res.pvalue, expected_pvalue, atol=1e-13)


Code review comments:
diff_hunk: @@ -1393,3 +1396,692 @@ def cramervonmises_2samp(x, y, method='auto'):
             p = max(0, 1. - _cdf_cvm_inf(tn))
 
     return CramerVonMisesResult(statistic=t, pvalue=p)
+
+
+attributes = ('statistic', 'pvalue', 'null_distribution')
+PermutationTestResult = make_dataclass('PermutationTestResult', attributes)
+
+
+# Functions like these are under review in gh-13312; nice to have them here
+# Remove when they are added to scipy/stats/stats.py or a utility file
+def _broadcast_arrays(arrays, axis=None):
+    """
+    Broadcast shapes of arrays, ignoring incompatibility of specified axes
+    """
+    new_shapes = _broadcast_array_shapes(arrays, axis=axis)
+    if axis is None:
+        new_shapes = [new_shapes]*len(arrays)
+    return [np.broadcast_to(array, new_shape)
+            for array, new_shape in zip(arrays, new_shapes)]
+
+
+def _broadcast_array_shapes(arrays, axis=None):
+    """
+    Broadcast shapes of arrays, ignoring incompatibility of specified axes
+    """
+    shapes = [np.asarray(arr).shape for arr in arrays]
+    return _broadcast_shapes(shapes, axis)
+
+
+def _broadcast_shapes(shapes, axis=None):
+    """
+    Broadcast shapes, ignoring incompatibility of specified axes
+    """
+    # First, ensure all shapes have same number of dimensions by prepending 1s.
+    n_dims = max([len(shape) for shape in shapes])
+    new_shapes = np.ones((len(shapes), n_dims), dtype=int)
+    for row, shape in zip(new_shapes, shapes):
+        row[len(row)-len(shape):] = shape  # can't use negative indices (-0:)
+
+    # Remove the shape elements of the axes to be ignored, but remember them.
+    if axis is not None:
+        axis = np.atleast_1d(axis)
+        axis[axis < 0] = n_dims + axis[axis < 0]
+        removed_shapes = new_shapes[:, axis]
+        new_shapes = np.delete(new_shapes, axis, axis=1)
+
+    # If arrays are broadcastable, shape elements that are 1 may be replaced
+    # with a corresponding non-1 shape element. Assuming arrays are
+    # broadcastable, that final shape element can be found with:
+    new_shape = np.max(new_shapes, axis=0)
+    # except in case of an empty array:
+    new_shape *= new_shapes.all(axis=0)
+
+    # Among all arrays, there can only be one unique non-1 shape element.
+    # Therefore, if any non-1 shape element does not match what we found
+    # above, the arrays must not be broadcastable after all.
+    if np.any(~((new_shapes == 1) | (new_shapes == new_shape))):
+        raise ValueError("Array shapes are incompatible for broadcasting.")
+
+    if axis is not None:
+        # Add back the shape elements that were ignored
+        new_axis = axis - np.arange(len(axis))
+        new_shapes = [tuple(np.insert(new_shape, new_axis, removed_shape))
+                      for removed_shape in removed_shapes]
+        return new_shapes
+    else:
+        return tuple(new_shape)
+
+
+def _all_partitions_concatenated(ns):
+    """
+    Generate all partitions of indices of groups of given sizes, concatenated
+    """
+    def all_partitions(z, n):
+        for c in combinations(z, n):
+            x0 = set(c)
+            x1 = z - x0
+            yield [x0, x1]
+
+    def all_partitions_n(z, ns):
+        if len(ns) == 0:
+            yield [z]
+            return
+        for c in all_partitions(z, ns[0]):
+            for d in all_partitions_n(c[1], ns[1:]):
+                yield c[0:1] + d
+
+    z = set(range(np.sum(ns)))
+    for partitioning in all_partitions_n(z, ns[:]):
+        x = np.concatenate([list(partition)
+                            for partition in partitioning]).astype(int)
+        yield x
+
+
+def _batch_generator(iterable, batch):
+    """A generator that yields batches of elements from an iterable"""
+    iterator = iter(iterable)
+    if batch <= 0:
+        raise ValueError("`batch` must be positive.")
+    while True:
+        z = []
+        try:
+            # get elements from iterator `batch` at a time
+            for i in range(batch):
+                z.append(next(iterator))
+            yield z
+        except StopIteration:
+            # when there are no more elements, yield the final batch and stop
+            if z:
+                yield z
+            break
+
+
+def _calculate_null_both(data, statistic, n_permutations, batch,
+                         random_state=None):
+    """
+    Calculate null distribution for independent sample tests.
+    """
+    n_samples = len(data)
+
+    # compute number of permutations
+    # (distinct partitions of data into samples of these sizes)
+    n_obs_i = [sample.shape[-1] for sample in data]  # observations per sample
+    n_obs_ic = np.cumsum(n_obs_i)
+    n_obs = n_obs_ic[-1]  # total number of observations
+    n_max = np.prod([comb(n_obs_ic[i], n_obs_ic[i-1])
+                     for i in range(n_samples-1, 0, -1)])
+
+    _memory_check(data, min(n_permutations, n_max))
+
+    # perm_generator is an iterator that produces permutations of indices
+    # from 0 to n_obs. We'll concatenate the samples, use these indices to
+    # permute the data, then split the samples apart again.
+    if n_permutations >= n_max:
+        n_permutations = n_max
+        perm_generator = _all_partitions_concatenated(n_obs_i)
+    else:
+        # Neither RandomState.permutation nor Generator.permutation
+        # can permute axis-slices independently. If this feature is
+        # added in the future, batches of the desired size should be
+        # generated in a single call.
+        perm_generator = (random_state.permutation(n_obs)
+                          for i in range(n_permutations))
+
+    batch = batch or int(n_permutations)
+    null_distribution = []
+
+    # First, concatenate all the samples. In batches, permute samples with
+    # indices produced by the `perm_generator`, split them into new samples of
+    # the original sizes, compute the statistic for each batch, and add these
+    # statistic values to the null distribution.
+    data = np.concatenate(data, axis=-1)
+    for indices in _batch_generator(perm_generator, batch=batch):
+        indices = np.array(indices)
+
+        # `indices` is 2D: each row is a permutation of the indices.
+        # We use it to index `data` along its last axis, which corresponds
+        # with observations.
+        # After indexing, the second to last axis of `data_batch` corresponds
+        # with permutations, and the last axis corresponds with observations.
+        data_batch = data[..., indices]
+
+        # Move the permutation axis to the front: we'll concatenate a list
+        # of batched statistic values along this zeroth axis to form the
+        # null distribution.
+        data_batch = np.moveaxis(data_batch, -2, 0)
+        data_batch = np.split(data_batch, n_obs_ic[:-1], axis=-1)
+        null_distribution.append(statistic(*data_batch, axis=-1))
+    null_distribution = np.concatenate(null_distribution, axis=0)
+
+    return null_distribution, n_permutations
+
+
+def _calculate_null_pairings(data, statistic, n_permutations, batch,
+                             random_state=None):
+    """
+    Calculate null distribution for association tests.
+    """
+    n_samples = len(data)
+
+    # compute number of permutations (factorial(n) permutations of each sample)
+    n_obs_sample = data[0].shape[-1]  # observations per sample; same for each
+    n_max = factorial(n_obs_sample)**n_samples
+
+    _memory_check(data, min(n_permutations, n_max))
+
+    # `perm_generator` is an iterator that produces a list of permutations of
+    # indices from 0 to n_obs_sample, one for each sample.
+    if n_permutations >= n_max:
+        n_permutations = n_max
+        # cartesian product of the sets of all permutations of indices
+        perm_generator = product(*(permutations(range(n_obs_sample))
+                                   for i in range(n_samples)))
+    else:
+        # Separate random permutations of indices for each sample.
+        # Again, it would be nice if RandomState/Generator.permutation
+        # could permute each axis-slice separately.
+        perm_generator = ([random_state.permutation(n_obs_sample)
+                           for i in range(n_samples)]
+                          for j in range(n_permutations))
+
+    batch = batch or int(n_permutations)
+    null_distribution = []
+
+    for indices in _batch_generator(perm_generator, batch=batch):
+        indices = np.array(indices)
+
+        # `indices` is 3D: the zeroth axis is for permutations, the next is
+        # for samples, and the last is for observations. Swap the first two
+        # to make the zeroth axis correspond with samples, as it does for
+        # `data`.
+        indices = np.swapaxes(indices, 0, 1)
+
+        # When we're done, `data_batch` will be a list of length `n_samples`.
+        # Each element will be a batch of random permutations of one sample.
+        # The zeroth axis of each batch will correspond with permutations,
+        # and the last will correspond with observations. (This makes it
+        # easy to pass into `statistic`.)
+        data_batch = [None]*n_samples
+        for i in range(n_samples):
+            data_batch[i] = data[i][..., indices[i]]
+            data_batch[i] = np.moveaxis(data_batch[i], -2, 0)
+
+        null_distribution.append(statistic(*data_batch, axis=-1))
+    null_distribution = np.concatenate(null_distribution, axis=0)
+
+    return null_distribution, n_permutations
+
+
+def _calculate_null_samples(data, statistic, n_permutations, batch,
+                            random_state=None):
+    """
+    Calculate null distribution for paired-sample tests.
+    """
+    n_samples = len(data)
+
+    # By convention, the meaning of the "samples" permutations type for
+    # data with only one sample is to flip the sign of the observations.
+    # Achieve this by adding a second sample - the negative of the original.
+    if n_samples == 1:
+        data = [data[0], -data[0]]
+
+    # The "samples" permutation strategy is the same as the "pairings"
+    # strategy except the roles of samples and observations are flipped.
+    # So swap these axes, then we'll use the function for the "pairings"
+    # strategy to do all the work!
+    data = np.swapaxes(data, 0, -1)
+
+    # (Of course, the user's statistic doesn't know what we've done here,
+    # so we need to pass it what it's expecting.)
+    def statistic_wrapped(*data, axis):
+        data = np.swapaxes(data, 0, -1)
+        if n_samples == 1:
+            data = data[0:1]
+        return statistic(*data, axis=axis)
+
+    return _calculate_null_pairings(data, statistic_wrapped, n_permutations,
+                                    batch, random_state)
+
+
+def _permutation_test_iv(data, statistic, permutation_type, vectorized,
+                         permutations, batch, alternative, axis, random_state):
+    """Input validation for `permutation_test`."""
+
+    axis_int = int(axis)
+    if axis != axis_int:
+        raise ValueError("`axis` must be an integer.")
+
+    permutation_types = {'samples', 'pairings', 'both'}
+    permutation_type = permutation_type.lower()
+    if permutation_type not in permutation_types:
+        raise ValueError(f"`permutation_type` must be in {permutation_types}.")
+
+    if vectorized not in {True, False}:
+        raise ValueError("`vectorized` must be `True` or `False`.")
+
+    message = "`data` must be a tuple containing at least two samples"
+    try:
+        if len(data) < 2 and permutation_type == "both":
+            raise ValueError(message)
+    except TypeError:
+        raise TypeError(message)
+
+    data = _broadcast_arrays(data, axis)
+    data_iv = []
+    for sample in data:
+        sample = np.atleast_1d(sample)
+        if sample.shape[axis] <= 1:
+            raise ValueError("each sample in `data` must contain two or more "
+                             "observations along `axis`.")
+        sample = np.moveaxis(sample, axis_int, -1)
+        data_iv.append(sample)
+
+    permutations_int = (int(permutations) if not np.isinf(permutations)
+                        else np.inf)
+    if permutations != permutations_int or permutations_int <= 0:
+        raise ValueError("`permutations` must be a positive integer.")
+
+    if batch is None:
+        batch_iv = batch
+    else:
+        batch_iv = int(batch)
+        if batch != batch_iv or batch_iv <= 0:
+            raise ValueError("`batch` must be a positive integer or None.")
+
+    alternatives = {'two-sided', 'greater', 'less'}
+    alternative = alternative.lower()
+    if alternative not in alternatives:
+        raise ValueError(f"`alternative` must be in {alternatives}")
+
+    random_state = check_random_state(random_state)
+
+    return (data_iv, statistic, permutation_type, vectorized, permutations_int,
+            batch_iv, alternative, axis_int, random_state)
+
+
+# This is mostly to protect me while developing. Eventually we want to batch
+# the computations, so this will be less of an issue, but if it sounds like a
+# good feature, let me know.
+def _memory_check(data, permutations):
+    needed = sum([sample.size for sample in data]) * permutations * 8
+    available = 1e9  # psutil.virtual_memory()[1] # can't use in SciPy
+    if needed > 0.8 * available:
+        raise MemoryError("Not enough memory available. "
+                          "Consider reducing the batch size.")
+
+
+def permutation_test(data, statistic, *, permutation_type='both',
+                     vectorized=False, permutations=np.inf, batch=None,
+                     alternative="two-sided", axis=0, random_state=None):
+    r"""
+    Performs a permutation test of a given statistic on provided data.
+
+    For independent sample statistics, the null hypothesis is that the data are
+    randomly sampled from the same distribution.
+    For paired sample statistics, two null hypothesis can be tested:
+    that the data are paired at random or that the data are assigned to samples
+    at random.
+
+    Parameters
+    ----------
+    data : iterable of array-like
+        Contains the samples, each of which is an array of observations.
+        Dimensions of sample arrays must be compatible for broadcasting except
+        along `axis`.
+    statistic : callable
+        Statistic for which the p-value of the hypothesis test is to be
+        calculated. `statistic` must be a callable that accepts samples
+        as separate arguments (e.g. ``statistic(*data)``) and returns the
+        resulting statistic.
+        If `vectorized` is set ``True``, `statistic` must also accept a keyword
+        argument `axis` and be vectorized to compute the statistic along the
+        provided `axis` of the sample arrays.
+    permutation_type : {'both', 'samples', 'pairings'}, optional
+        The type of permutations to be performed, in accordance with the
+        null hypothesis. The first two permutation types are for paired sample
+        statistics, in which all samples contain the same number of
+        observations and observations with corresponding indices along `axis`
+        are considered to be paired; the third is for independent sample
+        statistics.
+
+        - ``'samples'`` : observations are assigned to different samples
+          but remain paired with the same observations from other samples.
+          This permutation type is appropriate for paired sample hypothesis
+          tests such as the Wilcoxon signed-rank test and the paired t-test.
+        - ``'pairings'`` : observations are paired with different observations,
+          but they remain within the same sample. This permutation type is
+          appropriate for association/correlation tests with statistics such
+          as Spearman's :math:`\rho`, Kendall's :math:`\tau`, and Pearson's
+          :math:`r`.
+        - ``'both'`` (default) : observations are assigned to different samples
+          without preserving pairs. Samples may contain different numbers of
+          observations. This permutation type is appropriate for
+          independent sample hypothesis tests such as the Mann-Whitney
+          :math:`U` test and the independent sample t-test.
+
+          Please see the Notes section below for more detailed descriptions
+          of the permutation types.
+
+    vectorized : bool, default: ``False``
+        By default, `statistic` is assumed to calculate the statistic only for
+        1D arrays contained in `data`. If `vectorized` is set ``True``,
+        `statistic` must also accept a keyword argument `axis` and be
+        vectorized to compute the statistic along the provided `axis` of the ND
+        arrays in `data`. Use of a vectorized statistic can reduce computation
+        time.
+    permutations : int, default: ``np.inf``
path: scipy/stats/_hypotests.py
text: @author:
Something to consider before merging is renaming the `permutations` parameter to `n_resamples` to keep this consistent with `bootstrap` (and gh-14576). It would be different from the name of the equivalent parameter in `ttest_ind`, but I think it's more important to keep the names consistent among these "resampling methods". The plan (after this PR merges) would be to rename `_bootstrap.py` to `_resampling.py` and move this over there.

@user1:
any suggestion how to proceed? somehow `permutations` or `n_permutations` seems better than `n_resamples` since I would not directly associate the latter with permutations

@author:
Well, permutation testing is one of several ["resampling methods"](https://en.wikipedia.org/wiki/Resampling_(statistics)#Permutation_tests) in statistics, so I think `n_resamples` is definitely appropriate, and it would be more consistent with `bootstrap` and the proposed [`monte_carlo_test`](https://github.com/scipy/scipy/pull/14576) (other 
 resampling methods). You might begin to associate the two if the parameter were renamed; it would be clearly documented to mean the number of "permutations" (which, as we've discussed, is not a great name either) ; ) I'd suggest we either leave it as-is or change it to `n_resamples`. Your call!

@user1:
Ok, true, should we change it to `n_resamples` then?

@author:
Yes, I'll do that and change "both" to "independent" today!

@author:
Done.