Ctrl+K
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)[原始碼]#

對提供的資料執行給定統計量的置換檢定。

對於獨立樣本統計量，虛無假設是資料是從相同分佈中隨機抽樣而來。對於成對樣本統計量，可以檢定兩個虛無假設：資料是隨機成對的，或是資料是隨機分配到樣本中的。

參數:
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’s \(\rho\)、Kendall’s \(\tau\) 和 Pearson’s \(r\)。

  • 'independent' (預設) : 觀察值被分配到不同的樣本。樣本可能包含不同數量的觀察值。此置換類型適用於獨立樣本假設檢定，例如 Mann-Whitney \(U\) 檢定和獨立樣本 t 檢定。

    請參閱下方的「注意事項」章節，以取得置換類型的更詳細說明。

vectorizedbool, optional

如果 vectorized 設定為 False，則 statistic 不會傳遞關鍵字引數 axis，並且預期僅針對 1D 樣本計算統計量。如果 True，則 statistic 將會傳遞關鍵字引數 axis，並且預期在傳遞 ND 樣本陣列時，沿著 axis 計算統計量。如果 None (預設)，如果 axis 是 statistic 的參數，則 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] 中建議的不偏估計值。也就是說，在計算隨機虛無分佈中與檢定統計量的觀察值一樣極端比例時，分子和分母中的值都會加一。對此調整的一種解釋是，檢定統計量的觀察值始終包含在隨機虛無分佈的元素中。用於雙尾 p 值的慣例並非通用；如果偏好不同的定義，則傳回觀察到的檢定統計量和虛無分佈。

axisint, 預設值：0

要計算統計量的（廣播）樣本軸。如果樣本具有不同數量的維度，則在考慮 axis 之前，會將 singleton 維度預先加到維度較少的樣本中。

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

  • 如果 random_state 是 int，則會使用新的 RandomState 實例，並以 random_state 作為種子。

  • 如果 random_state 已經是 Generator 或 RandomState 實例，則會使用該實例。

在 1.15.0 版本中變更：作為從使用 numpy.random.RandomState 過渡到 numpy.random.Generator 的 SPEC-0007 轉換的一部分，此關鍵字已從 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 會以每個不同的方式正好置換一次。因此，針對 a 和 b 之間樣本的每個唯一配對，正好計算一次 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。令 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 會透過考慮虛無分佈中「接近」（在非精確 dtype 的浮點 epsilon 的 100 倍相對容差內）檢定統計量的觀察值的元素，將其視為等於檢定統計量的觀察值，來防止這種情況。但是，建議使用者檢查虛無分佈，以評估此比較方法是否合適，如果不合適，則手動計算 p 值。請參閱以下範例。

參考文獻

[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 會將虛無分佈中在 max(1e-14, abs(r)*1e-14) 範圍內的檢定統計量 r 的觀察值元素視為等於 r。

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

在大多數實際情況下，這種比較方法預期是準確的，但建議使用者透過檢查虛無分佈中接近統計量觀察值的元素來評估這一點。此外，請考慮使用可以使用精確算術計算的統計量（例如整數統計量）。

在本頁