Skip to content
Explore Automated Analysis

How to run a permutation test for a difference in means in SciPy (with Code)

Editorial · Published

The Short Answer (TL;DR)

Pass a difference-of-means function to scipy.stats.permutation_test with permutation_type="independent" for two exchangeable independent samples. Use an exact enumeration only when the number of group allocations is small, and state the direction and two-sided p-value convention.

import numpy as np; from scipy import stats
x, y = np.array([3.,5.,4.,6.]), np.array([1.,2.,4.,3.,2.])
r = stats.permutation_test((x,y), lambda a,b: a.mean()-b.mean(), permutation_type="independent", vectorized=False, n_resamples=np.inf, alternative="two-sided"); print(r.statistic, r.pvalue)

Detailed Implementation

This example uses synthetic data and does not need a network connection. Install the recorded package versions in a separate Python environment before running it.

python -m pip install numpy==2.5.4 scipy==1.18.1
from pathlib import Path
from math import comb
import numpy as np
from scipy import stats

# Independent synthetic observations in the same measurement units.
x = np.array([3., 5., 4., 6.])
y = np.array([1., 2., 4., 3., 2.])
if any(a.ndim != 1 or a.size < 2 or not np.isfinite(a).all() for a in (x,y)):
    raise ValueError("Use two finite one-dimensional samples with at least two values each.")
n_allocations = comb(x.size + y.size, x.size)
if n_allocations > 100_000:
    raise ValueError("Exact enumeration is too large for this example; design a seeded randomized test.")
def difference(a, b):
    return a.mean() - b.mean()
result = stats.permutation_test(
    (x,y), difference, permutation_type="independent", vectorized=False,
    n_resamples=np.inf, batch=126, alternative="two-sided",
)
report = (
    f"n_x={x.size}; n_y={y.size}; allocations={n_allocations}\n"
    f"mean_x-minus-mean_y={result.statistic:.6f} measurement units\n"
    f"exact two-sided p={result.pvalue:.9f}\n"
)
Path("permutation-means.txt").write_text(report, encoding="utf-8")
print(report)

The example reports 126 allocations, a difference of 2.100000 measurement units, and exact two-sided p=0.079365079.

What does the exact comparison enumerate?

There are 126 ways to assign four of the nine labelled observations to x and the remaining five to y. The observed difference is 2.1 measurement units. Equal numerical observations still occupy separate labelled positions, so repeated statistic values retain their multiplicity; deduplicating them would change their probabilities.

Which null hypothesis makes the labels exchangeable?

This independent-sample test assumes the observations can exchange group labels under a common-distribution null. Independence within and between groups matters. Equality of means alone does not establish exchangeability when the distributions or variances differ; an unstudentized mean-difference permutation test is not a universal unequal-variance mean test. A matched-pair design needs a different permutation scheme.

How are the two tails counted?

For the exact enumeration, SciPy doubles the smaller inclusive tail probability and caps the result at one. Here five of 126 allocations are at least as large as the observed statistic, giving 10/126. The other tail contains 125 allocations. This convention need not equal counting absolute values in an asymmetric null distribution.

What changes for a larger experiment?

The explicit guard prevents this small example from attempting a combinatorial workload. For a planned randomized approximation, choose a finite resample count and a recorded random generator seed, and report Monte Carlo uncertainty. SciPy adds one to the extreme count and denominator for randomized tests; the exact 126-allocation result here uses no such adjustment. Neither a small p-value nor the synthetic effect establishes causation.

Mathematical Context

The statistic is T = mean(x) minus mean(y), with the same units as the observations. Under exchangeability, each of the M = choose(n_x+n_y,n_x) labelled allocations is equally likely. Let L and G count null statistics at most and at least T_obs, including ties; the stated convention is p = min(1, 2 min(L/M, G/M)). In this example M=126, L=125 and G=5, so p=0.079365079 and the data do not reject the common-distribution null at a prespecified 0.05 level. A non-rejection does not prove the groups have the same distribution.

\[ p=\min\left(1,2\min\left(\frac{L}{M},\frac{G}{M}\right)\right) \]

Reproducibility

Executed with Python 3.12.14; numpy 2.5.4, scipy 1.18.1. Automated checks confirm the stated output properties.

Primary documentation

Related guides