scipy.stats.

permutation_test#

scipy.stats.permutation_test(data, statistic, *, permutation_type='independent', vectorized=None, n_resamples=9999, batch=None, alternative='two-sided', axis=0, rng=None, random_state=None)[源代码]#

对提供的数据执行给定统计量的排列检验（Permutation test）。

对于独立样本统计量，零假设是数据是从同一分布中随机抽样的。对于配对样本统计量，可以测试两个零假设：数据是随机配对的，或者数据是随机分配给样本的。

参数:
dataarray-like 的可迭代对象

包含样本，每个样本都是一个观测值数组。样本数组的维度必须兼容广播，但在 axis 轴上除外。

statisticcallable

需要计算假设检验 p 值的统计量。statistic 必须是一个可调用对象，接受样本作为独立参数（例如 statistic(*data)）并返回结果统计量。如果 vectorized 设置为 True，则 statistic 还必须接受关键字参数 axis，并支持向量化以沿样本数组提供的 axis 计算统计量。

permutation_type{‘independent’, ‘samples’, ‘pairings’}, 可选

根据零假设执行的排列类型。前两种排列类型用于配对样本统计量，其中所有样本包含相同数量的观测值，并且沿 axis 具有相应索引的观测值被视为配对；第三种用于独立样本统计量。

  • 'samples'：观测值被分配到不同的样本中，但仍与来自其他样本的相同观测值保持配对。这种排列类型适用于配对样本假设检验，如 Wilcoxon 符号秩检验和配对 t 检验。

  • 'pairings'：观测值与不同的观测值配对，但它们保留在同一个样本中。这种排列类型适用于关联/相关性测试，使用的统计量包括 Spearman’s \(\rho\)、Kendall’s \(\tau\) 和 Pearson’s \(r\)。

  • 'independent'（默认值）：观测值被分配到不同的样本中。样本可能包含不同数量的观测值。这种排列类型适用于独立样本假设检验，如 Mann-Whitney \(U\) 检验和独立样本 t 检验。

    请参阅下文的“注意”部分，了解排列类型的更详细说明。

vectorized布尔值，可选

如果 vectorized 设置为 False，statistic 将不会被传递关键字参数 axis，且预期仅计算 1D 样本的统计量。如果为 True，statistic 将被传递关键字参数 axis，且在传递 ND 样本数组时预期沿 axis 计算统计量。如果为 None（默认值），若 axis 是 statistic 的参数，则 vectorized 将被设置为 True。使用向量化统计量通常会减少计算时间。

n_resamplesint 或 np.inf, 默认值: 9999

用于近似零分布的随机排列（重采样）次数。如果大于或等于不同排列的总数，则将计算精确的零分布。请注意，不同排列的数量随样本大小增长非常迅速，因此精确检验仅适用于非常小的数据集。

batchint, 可选

每次调用 statistic 时处理的排列数量。内存使用量为 O( batch * n )，其中 n 是所有样本的总大小，无论 vectorized 的值如何。默认值为 None，在这种情况下 batch 即为排列的总数。

alternative{‘two-sided’, ‘less’, ‘greater’}, 可选

计算 p 值的备择假设。对于每种备择方案，精确检验的 p 值定义如下。

  • 'greater'：零分布中大于或等于观测统计量的百分比。

  • 'less'：零分布中小于或等于观测统计量的百分比。

  • 'two-sided'（默认值）：上述两个 p 值中较小者的两倍。

请注意，随机化检验的 p 值是根据 [2] 和 [3] 中建议的保守（高估）近似值计算的，而不是 [4] 中建议的无偏估计量。也就是说，在计算随机化零分布中与观测统计量同样极端的比例时，分子和分母的值都增加了 1。这种调整的一种解释是，观测到的测试统计量始终被视为随机化零分布的一个元素。双侧 p 值使用的约定并不通用；函数会返回观测统计量和零分布，以备需要不同的定义。

axisint, 默认值: 0

计算统计量的（广播后）样本轴。如果样本具有不同的维数，则在考虑 axis 之前，会将单一维度前置到维数较少的样本中。

rng{None, int, numpy.random.Generator}, 可选

如果通过关键字传递 rng，则 numpy.random.Generator 以外的类型将被传递给 numpy.random.default_rng 以实例化 Generator。如果 rng 已经是 Generator 实例，则使用提供的实例。指定 rng 以获得可重复的函数行为。

如果此参数通过位置传递，或通过关键字传递 random_state，则适用参数 random_state 的遗留行为

  • 如果 random_state 为 None（或 numpy.random），则使用 numpy.random.RandomState 单例。

  • 如果 random_state 是一个整数，则使用一个新的 RandomState 实例，并以 random_state 为种子。

  • 如果 random_state 已经是 Generator 或 RandomState 实例，则使用该实例。

1.15.0 版本中的变更：作为从使用 numpy.random.RandomState 向 numpy.random.Generator 过渡的 SPEC-007 的一部分，此关键字从 random_state 更改为 rng。在过渡期内，两个关键字都将继续有效，但一次只能指定一个。过渡期后，使用 random_state 关键字的函数调用将发出警告。random_state 和 rng 的行为如上所述，但在新代码中应仅使用 rng 关键字。

返回:
resPermutationTestResult

具有以下属性的对象

statisticfloat 或 ndarray

数据的观测测试统计量。

pvaluefloat 或 ndarray

给定备择假设下的 p 值。

null_distributionndarray

在零假设下生成的测试统计量的值。

附注

本函数支持的三种排列检验类型描述如下。

非配对统计量 (permutation_type='independent')

与此排列类型相关的零假设是所有观测值都采样自同一个潜在分布，并且它们被随机分配到其中一个样本中。

假设 data 包含两个样本；例如 a, b = data。当 1 < n_resamples < binom(n, k) 时，其中：

  • k 是 a 中的观测值数量，

  • n 是 a 和 b 中的观测值总数，且

  • binom(n, k) 是二项式系数（从 n 中取 k），

数据被池化（合并），随机分配到第一个或第二个样本中，并计算统计量。此过程重复执行 permutation 次，生成零假设下的统计量分布。将原始数据的统计量与此分布进行比较以确定 p 值。

当 n_resamples >= binom(n, k) 时，将执行精确检验：数据在样本之间以每种不同的方式进行一次划分，并形成精确的零分布。请注意，对于给定的样本间数据划分，仅考虑每个样本内部数据的一种排序/排列。对于不依赖于样本内数据顺序的统计量，这显著降低了计算成本而不影响零分布的形状（因为每个值的频率/计数都受到相同因子的影响）。

对于 a = [a1, a2, a3, a4] 和 b = [b1, b2, b3]，此类排列的一个示例是 x = [b3, a1, a2, b2] 和 y = [a4, b1, a3]。由于在精确检验中仅考虑每个样本内部数据的一种排序/排列，因此像 x = [b3, a1, b2, a2] 和 y = [a4, a3, b1] 这样的重采样将不被视为与上述示例不同。

permutation_type='independent' 不支持单样本统计量，但可以应用于具有两个以上样本的统计量。在这种情况下，如果 n 是每个样本中观测值数量的数组，则不同划分的数量为：

np.prod([binom(sum(n[i:]), sum(n[i+1:])) for i in range(len(n)-1)])

配对统计量，排列配对方式 (permutation_type='pairings')

与此排列类型相关的零假设是每个样本内的观测值都取自相同的潜在分布，并且与其他样本元素的配对是随机分配的。

假设 data 仅包含一个样本；例如 a, = data，我们希望考虑 a 的元素与第二个样本 b 的元素的所有可能配对。令 n 为 a 中的观测值数量，其也必须等于 b 中的观测值数量。

当 1 < n_resamples < factorial(n) 时，a 的元素被随机排列。用户提供的统计量接受一个数据参数，例如 a_perm，并考虑 a_perm 和 b 来计算统计量。此过程重复执行 permutation 次，生成零假设下的统计量分布。将原始数据的统计量与此分布进行比较以确定 p 值。

当 n_resamples >= factorial(n) 时，将执行精确检验：a 以每种不同的方式恰好排列一次。因此，statistic 将为 a 和 b 之间的每个唯一样本配对计算一次。

对于 a = [a1, a2, a3] 和 b = [b1, b2, b3]，此类排列的一个示例是 a_perm = [a3, a1, a2]，而 b 保持其原始顺序。

permutation_type='pairings' 支持包含任意数量样本的 data，每个样本必须包含相同数量的观测值。data 中提供的所有样本都是独立排列的。因此，如果 m 是样本数量，n 是每个样本中的观测值数量，则精确检验中的排列数量为：

factorial(n)**m

请注意，例如，如果一个双样本统计量本质上不依赖于提供观测值的顺序——而仅取决于观测值的配对——那么 data 中应仅提供两个样本中的一个。这显著降低了计算成本而不影响零分布的形状（因为每个值的频率/计数都受到相同因子的影响）。

配对统计量，排列样本 (permutation_type='samples')

与此排列类型相关的零假设是每对中的观测值都取自相同的潜在分布，并且它们被分配到的样本是随机的。

假设 data 包含两个样本；例如 a, b = data。令 n 为 a 中的观测值数量，其也必须等于 b 中的观测值数量。

当 1 < n_resamples < 2**n 时，a 和 b 的元素在样本之间随机交换（保持它们的配对关系）并计算统计量。此过程重复执行 permutation 次，生成零假设下的统计量分布。将原始数据的统计量与此分布进行比较以确定 p 值。

当 n_resamples >= 2**n 时，将执行精确检验：观测值以每种不同的方式分配到两个样本中（同时保持配对）各一次。

对于 a = [a1, a2, a3] 和 b = [b1, b2, b3]，此类排列的一个示例是 x = [b1, a2, b3] 和 y = [a1, b2, a3]。

permutation_type='samples' 支持包含任意数量样本的 data，每个样本必须包含相同数量的观测值。如果 data 包含多个样本，则 data 中的配对观测值在样本之间独立交换。因此，如果 m 是样本数量，n 是每个样本中的观测值数量，则精确检验中的排列数量为：

factorial(m)**n

几种配对样本统计检验，如 Wilcoxon 符号秩检验和配对样本 t 检验，可以仅考虑两个配对元素之间的差异。相应地，如果 data 仅包含一个样本，则零分布通过独立改变每个观测值的符号来形成。

警告

p 值是通过计算零分布中与观测到的统计量值同样极端或更极端的元素来计算的。由于使用有限精度算术，某些统计量函数在理论值完全相等时会返回数值上不同的值。在某些情况下，这可能会导致计算出的 p 值出现很大误差。permutation_test 通过将零分布中与观测到的测试统计量值“接近”（在非精确数据类型的浮点 epsilon 的 100 倍相对容差内）的元素视为等于观测到的测试统计量值来防止这种情况。但是，建议用户检查零分布以评估这种比较方法是否合适，如果不合适，则手动计算 p 值。请参见下面的示例。

数组 API 标准支持

permutation_test 除了支持 NumPy 之外，还对 Python Array API 标准兼容的后端提供实验性支持。请考虑通过设置环境变量 SCIPY_ARRAY_API=1 并提供 CuPy、PyTorch、JAX 或 Dask 数组作为数组参数来测试这些功能。支持以下后端和设备（或其他功能）的组合。

库

CPU

GPU

NumPy

✅

不适用

CuPy

不适用

✅

PyTorch

✅

✅

JAX

✅

✅

Dask

⛔

不适用

有关更多信息，请参阅 对数组 API 标准的支持。

参考文献

[1]

    1. Fisher. The Design of Experiments, 6th Ed (1951).

[2]

B. Phipson and G. K. Smyth. “Permutation P-values Should Never Be Zero: Calculating Exact P-values When Permutations Are Randomly Drawn.” Statistical Applications in Genetics and Molecular Biology 9.1 (2010).

[3]

M. D. Ernst. “Permutation Methods: A Basis for Exact Inference”. Statistical Science (2004).

[4]

B. Efron and R. J. Tibshirani. An Introduction to the Bootstrap (1993).

示例

假设我们希望测试两个样本是否来自同一分布。假设潜在分布对我们来说是未知的，并且在观察数据之前，我们假设第一个样本的均值会小于第二个样本的均值。我们决定使用样本均值之差作为测试统计量，并将 0.05 的 p 值视为具有统计学显著性。

为了提高效率，我们以向量化的方式编写定义测试统计量的函数：样本 x 和 y 可以是 ND 数组，并且将为沿 axis 的每个轴切片计算统计量。

>>> import numpy as np
>>> def statistic(x, y, axis):
...     return np.mean(x, axis=axis) - np.mean(y, axis=axis)

收集数据后，我们计算测试统计量的观测值。

>>> 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

事实上，测试统计量为负，这表明 x 背后分布的真实均值小于 y 背后分布的真实均值。为了确定如果两个样本来自同一分布，这种情况偶然发生的概率，我们执行排列检验。

>>> from scipy.stats import permutation_test
>>> # because our statistic is vectorized, we pass `vectorized=True`
>>> # `n_resamples=np.inf` indicates that an exact test is to be performed
>>> res = permutation_test((x, y), statistic, vectorized=True,
...                        n_resamples=np.inf, alternative='less')
>>> print(res.statistic)
-3.5411688580987266
>>> print(res.pvalue)
0.004329004329004329

在零假设下，获得小于或等于观测值的测试统计量的概率为 0.4329%。这小于我们选择的 5% 阈值，因此我们认为这是反对零假设、支持备择假设的显著证据。

由于上述样本的大小较小，permutation_test 可以执行精确检验。对于较大的样本，我们要求助于随机排列检验。

>>> x = norm.rvs(size=100, random_state=rng)
>>> y = norm.rvs(size=120, loc=0.2, random_state=rng)
>>> res = permutation_test((x, y), statistic, n_resamples=9999,
...                        vectorized=True, alternative='less',
...                        rng=rng)
>>> print(res.statistic)
-0.4230459671240913
>>> print(res.pvalue)
0.0015

在零假设下，获得小于或等于观测值的测试统计量的近似概率为 0.0225%。这再次小于我们选择的 5% 阈值，因此我们再次有显著证据拒绝零假设，转而支持备择假设。

对于大样本和大量排列，结果与相应的渐近检验（独立样本 t 检验）相当。

>>> from scipy.stats import ttest_ind
>>> res_asymptotic = ttest_ind(x, y, alternative='less')
>>> print(res_asymptotic.pvalue)
0.0014669545224902675

提供了测试统计量的排列分布，供进一步研究。

>>> 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")
>>> plt.show()
../../_images/scipy-stats-permutation_test-1_00_00.png

如果统计量由于机器精度有限而存在误差，则必须检查零分布。考虑以下情况：

>>> from scipy.stats import pearsonr
>>> x = [1, 2, 4, 3]
>>> y = [2, 4, 6, 8]
>>> def statistic(x, y, axis=-1):
...     return pearsonr(x, y, axis=axis).statistic
>>> res = permutation_test((x, y), statistic, vectorized=True,
...                        permutation_type='pairings',
...                        alternative='greater')
>>> r, pvalue, null = res.statistic, res.pvalue, res.null_distribution

在这种情况下，由于数值噪声，零分布的一些元素与相关系数 r 的观测值不同。我们手动检查零分布中与测试统计量观测值几乎相同的元素。

>>> r
0.7999999999999999
>>> unique = np.unique(null)
>>> unique
array([-1. , -1. , -0.8, -0.8, -0.8, -0.6, -0.4, -0.4, -0.2, -0.2, -0.2,
    0. ,  0.2,  0.2,  0.2,  0.4,  0.4,  0.6,  0.8,  0.8,  0.8,  1. ,
    1. ])  # may vary
>>> unique[np.isclose(r, unique)].tolist()
[0.7999999999999998, 0.7999999999999999, 0.8]  # may vary

如果 permutation_test 进行幼稚的比较，则值为 0.7999999999999998 的零分布元素将不会被视为与统计量观测值同样极端或更极端，因此计算出的 p 值将太小。

>>> incorrect_pvalue = np.count_nonzero(null >= r) / len(null)
>>> incorrect_pvalue
0.14583333333333334  # may vary

相反，permutation_test 将零分布中处于统计量观测值 r 的 max(1e-14, abs(r)*1e-14) 范围内的元素视为等于 r。

>>> correct_pvalue = np.count_nonzero(null >= r - 1e-14) / len(null)
>>> correct_pvalue
0.16666666666666666
>>> res.pvalue == correct_pvalue
True

这种比较方法在大多数实际情况下预计是准确的，但建议用户通过检查零分布中接近统计量观测值的元素来评估这一点。此外，考虑使用可以使用精确算术计算的统计量（例如整数统计量）。

本页内容