monte_carlo_test#
- scipy.stats.monte_carlo_test(data, rvs, statistic, *, vectorized=None, n_resamples=9999, batch=None, alternative='two-sided', axis=0)[原始碼]#
執行蒙地卡羅假設檢定。
data 包含一個樣本或一個或多個樣本的序列。rvs 指定了在虛無假設下 data 中樣本的分佈。statistic 對於給定 data 的值會與蒙地卡羅虛無分佈進行比較:使用 rvs 生成的每組 n_resamples 樣本的統計量值。這給出了 p 值,即在虛無假設下觀察到如此極端檢定統計量值的機率。
- 參數:
- dataarray-like 或 array-like 序列
觀測值的陣列或陣列序列。
- rvs可呼叫物件或可呼叫物件的元組
一個可呼叫物件或可呼叫物件的序列,用於在虛無假設下生成隨機變數。 rvs 的每個元素都必須是一個可呼叫物件,它接受關鍵字引數
size
(例如rvs(size=(m, n))
) 並返回該形狀的 N 維陣列樣本。如果 rvs 是一個序列,則 rvs 中可呼叫物件的數量必須與 data 中的樣本數量相符,即len(rvs) == len(data)
。如果 rvs 是一個單一可呼叫物件,則 data 被視為單一樣本。- statistic可呼叫物件
要計算假設檢定 p 值的統計量。statistic 必須是一個可呼叫物件,它接受一個樣本 (例如
statistic(sample)
) 或len(rvs)
個獨立樣本 (例如,如果 rvs 包含兩個可呼叫物件且 data 包含兩個樣本,則為statistic(samples1, sample2)
),並返回結果統計量。如果 vectorized 設定為True
,則 statistic 也必須接受關鍵字引數 axis,並且可以向量化以沿著 data 中樣本的提供的 axis 計算統計量。- vectorizedbool,選用
如果 vectorized 設定為
False
,則不會將關鍵字引數 axis 傳遞給 statistic,並且預期僅針對 1D 樣本計算統計量。如果True
,則會將關鍵字引數 axis 傳遞給 statistic,並且預期在傳遞 ND 樣本陣列時沿著 axis 計算統計量。如果None
(預設值),如果axis
是 statistic 的參數,則 vectorized 將設定為True
。使用向量化統計量通常可以減少計算時間。- n_resamplesint,預設值:9999
從每個 rvs 的可呼叫物件中抽取的樣本數。等效地,虛無假設下用作蒙地卡羅虛無分佈的統計量值數量。
- batchint,選用
每次呼叫 statistic 時要處理的蒙地卡羅樣本數。記憶體用量為 O( batch *
sample.size[axis]
)。預設值為None
,在這種情況下,batch 等於 n_resamples。- alternative{‘two-sided’, ‘less’, ‘greater’}
計算 p 值的對立假設。對於每個對立假設,p 值定義如下。
'greater'
:虛無分佈中大於或等於檢定統計量的觀測值的百分比。'less'
:虛無分佈中小於或等於檢定統計量的觀測值的百分比。'two-sided'
:上述 p 值中較小值的兩倍。
- axisint,預設值:0
data (或 data 中每個樣本) 中要計算統計量的軸。
- 返回:
- resMonteCarloTestResult
具有以下屬性的物件
- statisticfloat 或 ndarray
觀測到的 data 的檢定統計量。
- pvaluefloat 或 ndarray
給定對立假設的 p 值。
- null_distributionndarray
在虛無假設下生成的檢定統計量的值。
警告
p 值的計算方式是計算虛無分佈中與統計量的觀測值一樣極端或更極端的元素數量。由於使用了有限精確度算術,因此當理論值完全相等時,某些統計函數會返回數值上不同的值。在某些情況下,這可能會導致計算出的 p 值出現較大誤差。
monte_carlo_test
通過考慮虛無分佈中「接近」(在不精確資料類型的浮點 epsilon 的 100 倍相對容差內)檢定統計量的觀測值的元素,將它們視為等於檢定統計量的觀測值,以此來防止這種情況。但是,建議使用者檢查虛無分佈以評估這種比較方法是否合適,如果不合適,則手動計算 p 值。
參考文獻
[1]B. Phipson 和 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)。
範例
假設我們希望檢定小樣本是否是從常態分佈中抽取的。我們決定使用樣本的偏度作為檢定統計量,並且我們將認為 p 值為 0.05 在統計上是顯著的。
>>> import numpy as np >>> from scipy import stats >>> def statistic(x, axis): ... return stats.skew(x, axis)
在收集我們的資料後,我們計算檢定統計量的觀測值。
>>> rng = np.random.default_rng() >>> x = stats.skewnorm.rvs(a=1, size=50, random_state=rng) >>> statistic(x, axis=0) 0.12457412450240658
為了確定如果樣本是從常態分佈中抽取,則通過偶然性觀察到如此極端偏度的機率,我們可以執行蒙地卡羅假設檢定。該檢定將從常態分佈中隨機抽取許多樣本,計算每個樣本的偏度,並將我們的原始偏度與此分佈進行比較,以確定近似的 p 值。
>>> from scipy.stats import monte_carlo_test >>> # because our statistic is vectorized, we pass `vectorized=True` >>> rvs = lambda size: stats.norm.rvs(size=size, random_state=rng) >>> res = monte_carlo_test(x, rvs, statistic, vectorized=True) >>> print(res.statistic) 0.12457412450240658 >>> print(res.pvalue) 0.7012
在虛無假設下獲得小於或等於觀測值的檢定統計量的機率約為 70%。這大於我們選擇的 5% 閾值,因此我們不能認為這是反對虛無假設的顯著證據。
請注意,此 p 值基本上與
scipy.stats.skewtest
的 p 值相符,後者依賴於基於樣本偏度的檢定統計量的漸近分佈。>>> stats.skewtest(x).pvalue 0.6892046027110614
這種漸近近似對於小樣本量無效,但
monte_carlo_test
可以用於任何大小的樣本。>>> x = stats.skewnorm.rvs(a=1, size=7, random_state=rng) >>> # stats.skewtest(x) would produce an error due to small sample >>> res = monte_carlo_test(x, rvs, statistic, vectorized=True)
檢定統計量的蒙地卡羅分佈已提供,以供進一步研究。
>>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots() >>> ax.hist(res.null_distribution, bins=50) >>> ax.set_title("Monte Carlo distribution of test statistic") >>> ax.set_xlabel("Value of Statistic") >>> ax.set_ylabel("Frequency") >>> plt.show()