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)[source]#

对提供的数据上的给定统计量执行置换检验。

对于独立样本统计量,零假设是数据是从同一分布中随机抽样的。对于配对样本统计量,可以检验两个零假设:数据是随机配对的,或者数据是随机分配到样本中的。

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

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

statisticcallable

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

permutation_type{‘independent’, ‘samples’, ‘pairings’}, optional

要执行的置换的类型,与零假设一致。前两个置换类型用于配对样本统计量,其中所有样本包含相同数量的观测值,并且沿 axis 具有对应索引的观测值被认为是配对的;第三个用于独立样本统计量。

  • 'samples': 观测值被分配到不同的样本,但与其他样本的相同观测值保持配对。此置换类型适用于配对样本假设检验,例如 Wilcoxon 符号秩检验和配对 t 检验。

  • 'pairings': 观测值与不同的观测值配对,但它们保留在同一样本中。此置换类型适用于具有统计量的关联/相关性检验,例如 Spearman 的 \(\rho\)、Kendall 的 \(\tau\) 和 Pearson 的 \(r\)

  • 'independent' (默认): 观测值被分配到不同的样本。样本可能包含不同数量的观测值。此置换类型适用于独立样本假设检验,例如 Mann-Whitney \(U\) 检验和独立样本 t 检验。

    请参阅下面的“注释”部分,了解置换类型的更详细说明。

vectorizedbool, optional

如果 vectorized 设置为 False,则不会将关键字参数 axis 传递给 statistic,并且期望仅为 1D 样本计算统计量。如果 True,则会将关键字参数 axis 传递给 statistic,并且期望在传递 ND 样本数组时沿 axis 计算统计量。如果 None (默认),如果 axisstatistic 的参数,则 vectorized 将设置为 True。使用向量化统计量通常会减少计算时间。

n_resamplesint 或 np.inf,默认值: 9999

用于近似零分布的随机置换(重采样)的数量。如果大于或等于不同置换的数量,则将计算精确的零分布。请注意,不同置换的数量随着样本大小的增加而迅速增加,因此精确检验仅适用于非常小的数据集。

batchint, optional

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

alternative{‘two-sided’, ‘less’, ‘greater’}, optional

计算 p 值的备择假设。对于每个备择假设,p 值定义为精确检验如下。

  • 'greater': 大于或等于检验统计量的观测值的零分布的百分比。

  • 'less': 小于或等于检验统计量的观测值的零分布的百分比。

  • 'two-sided' (默认): 上面 p 值中较小的两倍。

请注意,随机检验的 p 值是根据 [2][3] 中建议的保守(高估)近似值计算的,而不是 [4] 中建议的无偏估计量。也就是说,当计算随机零分布中与检验统计量的观测值一样极端的比例时,分子和分母中的值都加 1。对此调整的一种解释是,检验统计量的观测值始终包含为随机零分布的一个元素。用于双侧 p 值的约定不是通用的;观测到的检验统计量和零分布将返回,以防需要不同的定义。

axisint, default: 0

在 (广播) 样本上计算统计量的轴。如果样本具有不同数量的维度,则在考虑 axis 之前,将单例维度添加到具有较少维度的样本。

rng{None, int, numpy.random.Generator}, optional

如果通过关键字传递 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 已经是 GeneratorRandomState 实例,则使用该实例。

在 1.15.0 版本中更改: 作为从使用 numpy.random.RandomState 转换为 numpy.random.GeneratorSPEC-007 过渡的一部分,此关键字已从 random_state 更改为 rng。在过渡期间,这两个关键字将继续工作,尽管一次只能指定一个关键字。在过渡期之后,使用 random_state 关键字的函数调用将发出警告。上面概述了 random_staterng 的行为,但新代码中应仅使用 rng 关键字。

返回:
resPermutationTestResult

具有属性的对象

statisticfloat 或 ndarray

数据的观测检验统计量。

pvaluefloat 或 ndarray

给定备择假设的 p 值。

null_distributionndarray

在零假设下生成的检验统计量的值。

注释

下面介绍了此函数支持的三种置换检验。

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

与此置换类型关联的零假设是所有观测值都从同一底层分布中抽样,并且它们已随机分配给其中一个样本。

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

  • ka 中的观测值数量,

  • nab 中的观测值总数,并且

  • 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 的元素的所有可能配对。令 na 中的观测值数量,该数量也必须等于 b 中的观测值数量。

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

n_resamples >= factorial(n) 时,将执行精确检验:a 以每个不同的方式仅置换一次。因此,对于 ab 之间样本的每个唯一配对,statistic 仅计算一次。

对于 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。令 na 中的观测值数量,该数量也必须等于 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 通过考虑零分布中“接近”(在不精确的 dtype 的浮点 epsilon 的 100 倍的相对容差范围内)于检验统计量的观测值的元素与检验统计量的观测值相等来防止这种情况。但是,建议用户检查零分布以评估此比较方法是否合适,如果不是,则手动计算 p 值。请参见下面的示例。

参考文献

[1]
    1. Fisher. 实验设计,第 6 版 (1951)。

[2]

B. Phipson 和 G. K. Smyth。“排列 P 值绝不应为零:当排列是随机绘制时计算精确 P 值。” 遗传学和分子生物学统计应用 9.1 (2010)。

[3]

M. D. Ernst。“排列方法:精确推断的基础”。 统计科学 (2004)。

[4]

B. Efron 和 R. J. Tibshirani。 Bootstrap 简介 (1993)。

示例

假设我们希望测试两个样本是否来自相同的分布。 假设我们不知道底层分布,并且在观察数据之前,我们假设第一个样本的平均值将小于第二个样本的平均值。 我们决定使用样本均值之间的差异作为检验统计量,并且我们将认为 0.05 的 p 值为统计显着。

为了提高效率,我们以矢量化的方式编写定义检验统计量的函数:样本 xy 可以是 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

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