scipy.stats.

power#

scipy.stats.power(test, rvs, n_observations, *, significance=0.01, vectorized=None, n_resamples=10000, batch=None, kwargs=None)[源代码]#

在替代假设下模拟假设检验的能力。

参数:
test可调用

对其计算功效的假设检验。 test 必须是可以接受样本的 callable(例如 test(sample))或 len(rvs) 个独立样本(例如 test(samples1, sample2) 其中 rvs 包含两个 callable, n_observations 包含两个值)并返回检验的 p 值。如果将 vectorized 设为 True,则 test 也必须可以接收关键字参数 axis,并矢量化以沿样本中提供的 axis 执行检验。scipy.stats 中带有 axis 参数(该参数返回带有 pvalue 属性的对象)的任何 callable 也可以接受。

rvscallable 或 callable 元组

一个 callable 或可生成备择假设下随机变量的 callable 序列。 rvs 的每个元素都必须接受关键字参数 size(例如 rvs(size=(m, n))),并返回相应形状的 N 维数组。如果 rvs 是一个序列,那么 rvs 中的 callable 数量必须与 n_observations 的元素数量相匹配,即 len(rvs) == len(n_observations)。如果 rvs 是单个 callable,则 n_observations 被当作一个元素处理。

n_observations整数元组或整数数组元组

如果是一个整数序列,则每个整数都是传递给 test 的样本的大小。如果是一个整数数组序列,则会针对每组对应的样本大小模拟功效。请参见示例。

significance浮点数或浮点数数组_类似物,默认值:0.01

显著性阈值;即,将假设检验结果视为反对无效假设的证据之下的p值。等价地即空假设下的 I 类错误的可接受比率。如果是一个数组,则会针对每个显著性阈值模拟功效。

kwargs可选,字典

传递给 rvs 和/或 test 可调用对象的关键字参数。使用自省来确定哪些关键字参数可以传递给每个可调用对象。每个关键字对应的值必须是一个数组。数组必须彼此广播,并与 n_observations 中的每个数组广播。针对每组对应的样本大小和参数模拟幂。参见示例。

vectorizedbool,可选

如果 vectorized 设置为 False,则不会将关键字参数 axis 传递给 test,并且仅针对一维样本执行测试。如果 True,则关键字参数 axis 将传递给 test,并且在传递 N 维样本数组时,预期沿 axis 执行测试。如果 None(默认),则当 axistest 的参数时,vectorized 将设置为 True。通常,使用矢量化检验可以减少计算时间。

n_resamplesint,默认值:10000

rvs 的每个可调用对象中提取的样本数量。同样,为了近似幂,在备择假设下执行的测试数量。

batchint,可选

在每次调用 test 时处理的样本数量。内存使用量与 batch 与最大样本大小的乘积成正比。默认值为 None,在这种情况下,batch 等于 n_resamples

返回:
resPowerResult

具有以下属性的对象

powerfloat 或 ndarray

针对备择估计的幂。

pvaluesndarray

在备择假设下观察到的 p 值。

注释

幂值模拟如下

  • rvs 指定的备选方案下,使用 n_observations 指定的大小(或一组样本的大小)绘制许多随机样本(或一组样本)。

  • 针对每个样本(或一组样本),根据 test 计算 p 值。这些 p 值记录在结果对象 pvalues 属性中。

  • 计算小于 significance 显著性水平的 p 值所占的比例。这是记录在结果对象 power 属性的幂值。

假定 significance 是一个形状为 shape1 的数组,kwargsn_observations 的元素可以相互广播到形状 shape2,并且 test 输出一个形状为 shape3 的 p 值数组。然后结果对象 power 属性形状将为 shape1 + shape2 + shape3pvalues 属性形状将为 shape2 + shape3 + (n_resamples,)

示例

假设我们希望在以下情况下模拟独立样本 t 检验的幂值

  • 第一个样本有 10 个观察值,从均值为 0 的正态分布中提取得到。

  • 第二个样本有 12 个观察值,从均值为 1.0 的正态分布中提取得到。

  • 显著性 p 值的阈值为 0.05。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>>
>>> test = stats.ttest_ind
>>> n_observations = (10, 12)
>>> rvs1 = rng.normal
>>> rvs2 = lambda size: rng.normal(loc=1, size=size)
>>> rvs = (rvs1, rvs2)
>>> res = stats.power(test, rvs, n_observations, significance=0.05)
>>> res.power
0.6116

在样本量分别为 10 和 12 的情况下,显著性阈值为 0.05 的 t 检验的幂值在所选备选方案下约为 60%。我们可以通过传递样本大小数组来研究样本大小对幂值的影响。

>>> import matplotlib.pyplot as plt
>>> nobs_x = np.arange(5, 21)
>>> nobs_y = nobs_x
>>> n_observations = (nobs_x, nobs_y)
>>> res = stats.power(test, rvs, n_observations, significance=0.05)
>>> ax = plt.subplot()
>>> ax.plot(nobs_x, res.power)
>>> ax.set_xlabel('Sample Size')
>>> ax.set_ylabel('Simulated Power')
>>> ax.set_title('Simulated Power of `ttest_ind` with Equal Sample Sizes')
>>> plt.show()
../../_images/scipy-stats-power-1_00_00.png

还可以研究效应大小对幂值的影响。在本例中,效应大小是第二个样本的分布位置。

>>> n_observations = (10, 12)
>>> loc = np.linspace(0, 1, 20)
>>> rvs2 = lambda size, loc: rng.normal(loc=loc, size=size)
>>> rvs = (rvs1, rvs2)
>>> res = stats.power(test, rvs, n_observations, significance=0.05,
...                   kwargs={'loc': loc})
>>> ax = plt.subplot()
>>> ax.plot(loc, res.power)
>>> ax.set_xlabel('Effect Size')
>>> ax.set_ylabel('Simulated Power')
>>> ax.set_title('Simulated Power of `ttest_ind`, Varying Effect Size')
>>> plt.show()
../../_images/scipy-stats-power-1_01_00.png

我们还可以使用 power 来估计检验的 I 类错误率(也被模棱两可的术语“大小”指代),并评估它是否与标称水平相匹配。例如,jarque_bera 的原假设是样本取自与正态分布具有相同偏度和峰度的分布。为了估计 I 类错误率,我们可以将原假设视为真备择假设并计算能力。

>>> test = stats.jarque_bera
>>> n_observations = 10
>>> rvs = rng.normal
>>> significance = np.linspace(0.0001, 0.1, 1000)
>>> res = stats.power(test, rvs, n_observations, significance=significance)
>>> size = res.power

如下所示,该检验的 I 类错误率远远低于其文档中提到的,如此小的样本的标称水平。

>>> ax = plt.subplot()
>>> ax.plot(significance, size)
>>> ax.plot([0, 0.1], [0, 0.1], '--')
>>> ax.set_xlabel('nominal significance level')
>>> ax.set_ylabel('estimated test size (Type I error rate)')
>>> ax.set_title('Estimated test size vs nominal significance level')
>>> ax.set_aspect('equal', 'box')
>>> ax.legend(('`ttest_1samp`', 'ideal test'))
>>> plt.show()
../../_images/scipy-stats-power-1_02_00.png

正如人们对如此保守的检验所期望的那样,对于一些备择假设,能力相当低。例如,检验在备择假设下(样本取自拉普拉斯分布)的能力可能不会比 I 类错误率大多少。

>>> rvs = rng.laplace
>>> significance = np.linspace(0.0001, 0.1, 1000)
>>> res = stats.power(test, rvs, n_observations, significance=0.05)
>>> print(res.power)
0.0587

这不是 SciPy 的实现错误;它只是由于以下事实,即检验统计量的原分布是在样本量很大的假设(即接近无穷大)下得出的,并且这种渐近近似对于小样本来说是不准确的。在这种情况下,重新采样和蒙特卡罗方法(例如,permutation_testgoodness_of_fitmonte_carlo_test)可能更合适。