SciPy-1-12-中文文档-二十六-

SciPy 1.12 中文文档(二十六)

原文:docs.scipy.org/doc/scipy-1.12.0/index.html

scipy.stats.chisquare

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.chisquare.html#scipy.stats.chisquare

scipy.stats.chisquare(f_obs, f_exp=None, ddof=0, axis=0)

计算单向卡方检验。

卡方检验检验分类数据是否具有给定频率的零假设。

参数:

f_obsarray_like

每个类别中的观察频率。

f_exparray_like,可选

每个类别中的期望频率。默认情况下,假定类别是等可能的。

ddofint,可选

“Δ自由度”:用于 p 值的自由度调整。p 值使用具有k - 1 - ddof自由度的卡方分布计算,其中k是观察频率的数量。ddof的默认值为 0。

axisint 或 None,可选

广播结果的轴f_obsf_exp在其上应用测试。如果 axis 为 None,则将f_obs中的所有值视为单个数据集。默认为 0。

返回:

res:Power_divergenceResult

包含属性的对象:

statisticfloat 或 ndarray

卡方检验统计量。如果 axis 为 None 或f_obsf_exp为 1-D,则该值为浮点数。

pvaluefloat 或 ndarray

测试的 p 值。如果ddof和结果属性statistic是标量,则该值为浮点数。

参见

scipy.stats.power_divergence

scipy.stats.fisher_exact

2x2 列联表上的 Fisher 确切性检验。

scipy.stats.barnard_exact

无条件精确性检验。对小样本量的卡方检验的替代方法。

注释

当每个类别中的观察或期望频率太小时,此检验无效。一个典型的规则是所有观察和期望频率应至少为 5。根据[3],推荐总样本数大于 13,否则应使用精确测试(如巴纳德精确检验),因为它们不会过度拒绝。

此外,观察频率和期望频率的总和必须相同才能使测试有效;如果在相对容差为1e-8的情况下这些总和不一致,chisquare 将会引发错误。

默认的自由度,k-1,适用于在未估计分布参数的情况下。如果通过高效的最大似然估计了 p 个参数,则正确的自由度为 k-1-p。如果参数以不同方式估计,则自由度可以在 k-1-p 和 k-1 之间。然而,也有可能渐近分布不是卡方分布,此时这个检验就不合适了。

参考文献

[1]

Lowry, Richard。“推断统计学的概念与应用”。第八章。web.archive.org/web/20171022032306/http://vassarstats.net:80/textbook/ch8pt1.html

[2]

“卡方检验”,zh.wikipedia.org/wiki/卡方检验

[3]

Pearson, Karl。“关于假设,即在相关系统中,给定的偏差系统的准则是这样的,可以合理地假设其是由随机抽样产生的”,《哲学杂志》。第 5 系列。50 (1900),第 157-175 页。

[4]

Mannan, R. William 和 E. Charles. Meslow. “俄勒冈东北部管理和原始森林中的鸟类种群和植被特征”。《野生动物管理杂志》48,1219-1238,DOI:10.2307/3801783,1984 年。

示例

[4] 中,研究了俄勒冈州一片古老的原始森林中的鸟类觅食行为。在这片森林中,有 44% 的冠层体积是道格拉斯冷杉,24% 是黄松,29% 是大冷杉,3% 是西部落叶松。作者观察了几种鸟类的行为,其中之一是红胁䴓。他们对这种鸟类的觅食行为进行了 189 次观察,并记录了在观察中道格拉斯冷杉中有 43(“23%”),黄松中有 52(“28%”),大冷杉中有 54(“29%”),西部落叶松中有 40(“21%”)次观察。

使用卡方检验,我们可以测试零假设,即觅食事件的比例等于树冠层体积的比例。文章的作者认为 p 值小于 1%是显著的。

使用上述冠层体积和观察事件的比例,我们可以推断期望频率。

>>> import numpy as np
>>> f_exp = np.array([44, 24, 29, 3]) / 100 * 189 

观察到的觅食频率为:

>>> f_obs = np.array([43, 52, 54, 40]) 

现在我们可以将观察频率与期望频率进行比较。

>>> from scipy.stats import chisquare
>>> chisquare(f_obs=f_obs, f_exp=f_exp)
Power_divergenceResult(statistic=228.23515947653874, pvalue=3.3295585338846486e-49) 

p 值远低于选定的显著水平。因此,作者认为差异显著,并得出结论,觅食事件的相对比例与树冠层体积的相对比例不同。

以下是其他通用示例,用于演示如何使用其他参数。

当只给出 f_obs 时,假定期望频率是均匀的,并由观察频率的平均值给出。

>>> chisquare([16, 18, 16, 14, 12, 12])
Power_divergenceResult(statistic=2.0, pvalue=0.84914503608460956) 

使用 f_exp 可以提供期望频率。

>>> chisquare([16, 18, 16, 14, 12, 12], f_exp=[16, 16, 16, 16, 16, 8])
Power_divergenceResult(statistic=3.5, pvalue=0.62338762774958223) 

f_obs 是 2-D 时,默认情况下将测试应用于每一列。

>>> obs = np.array([[16, 18, 16, 14, 12, 12], [32, 24, 16, 28, 20, 24]]).T
>>> obs.shape
(6, 2)
>>> chisquare(obs)
Power_divergenceResult(statistic=array([2\.        , 6.66666667]), pvalue=array([0.84914504, 0.24663415])) 

通过设置 axis=None,可以将测试应用于数组中的所有数据,这相当于将测试应用于展平的数组。

>>> chisquare(obs, axis=None)
Power_divergenceResult(statistic=23.31034482758621, pvalue=0.015975692534127565)
>>> chisquare(obs.ravel())
Power_divergenceResult(statistic=23.310344827586206, pvalue=0.01597569253412758) 

ddof 是对默认自由度的更改。

>>> chisquare([16, 18, 16, 14, 12, 12], ddof=1)
Power_divergenceResult(statistic=2.0, pvalue=0.7357588823428847) 

通过使用 ddof 广播卡方统计量来计算 p 值。

>>> chisquare([16, 18, 16, 14, 12, 12], ddof=[0,1,2])
Power_divergenceResult(statistic=2.0, pvalue=array([0.84914504, 0.73575888, 0.5724067 ])) 

f_obsf_exp 也会进行广播。在下面的例子中,f_obs 的形状为 (6,),f_exp 的形状为 (2, 6),因此广播 f_obsf_exp 的结果形状为 (2, 6)。为了计算所需的卡方统计量,我们使用 axis=1

>>> chisquare([16, 18, 16, 14, 12, 12],
...           f_exp=[[16, 16, 16, 16, 16, 8], [8, 20, 20, 16, 12, 12]],
...           axis=1)
Power_divergenceResult(statistic=array([3.5 , 9.25]), pvalue=array([0.62338763, 0.09949846])) 

scipy.stats.power_divergence

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.power_divergence.html#scipy.stats.power_divergence

scipy.stats.power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None)

Cressie-Read 功效散度统计量和拟合优度检验。

该函数使用 Cressie-Read 功效散度统计量检验分类数据具有给定频率的零假设。

参数:

f_obs:类数组

每个类别中的观察频率。

f_exp:类数组,可选

每个类别中的期望频率。默认情况下,假定类别是等可能的。

ddof:整数,可选

“Delta 自由度”:调整 p 值的自由度。使用自由度为k - 1 - ddof的卡方分布计算 p 值,其中k为观察频率的数量。ddof的默认值为 0。

axis:整数或 None,可选

沿着其应用测试的f_obsf_exp的广播结果的轴。如果轴为 None,则所有f_obs中的值都视为单个数据集。默认为 0。

lambda_:浮点数或字符串,可选

Cressie-Read 功效散度统计量的功率。默认值为 1。为方便起见,lambda_可以分配以下字符串之一,此时将使用相应的数值:

  • "pearson"(值为 1)

    Pearson 的卡方统计量。在这种情况下,该函数等同于chisquare

  • "log-likelihood"(值为 0)

    对数似然比。也称为 G 检验[3]

  • "freeman-tukey"(值为-1/2)

    Freeman-Tukey 统计量。

  • "mod-log-likelihood"(值为-1)

    修改的对数似然比。

  • "neyman"(值为-2)

    Neyman 统计量。

  • "cressie-read"(值为 2/3)

    推荐的功率见[5]

返回:

res:Power_divergenceResult

包含以下属性的对象:

统计量:浮点数或数组

Cressie-Read 功效散度检验统计量。如果axis为 None 或f_obsf_exp为 1-D,则该值为浮点数。

p 值:浮点数或数组

测试的 p 值。如果ddof和返回值stat为标量,则该值为浮点数。

另请参见

chisquare

注意

当每个类别中的观察或期望频率过小时,该检验无效。通常规则是所有观察和期望频率都应至少为 5。

此外,测试有效时观察和期望频率的总和必须相同;如果不同意则power_divergence会引发错误,相对容差为1e-8

lambda_ 小于零时,统计量的公式涉及除以 f_obs,因此如果 f_obs 中的任何值为零,则可能生成警告或错误。

类似地,如果在 lambda_ >= 0 时 f_exp 中的任何值为零,可能会生成警告或错误。

默认的自由度 k-1 适用于分布参数未估计的情况。如果通过高效的最大似然估计估计了 p 个参数,则正确的自由度为 k-1-p。如果以不同的方式估计参数,则自由度可以在 k-1-p 和 k-1 之间。然而,也有可能渐近分布不是卡方分布,在这种情况下,此检验不适用。

此函数处理屏蔽数组。如果 f_obsf_exp 的元素被屏蔽,则忽略该位置的数据,并且不计入数据集的大小。

新版本 0.13.0 中引入。

参考资料

[1]

Lowry, Richard。“推断统计学的概念与应用”。第八章。web.archive.org/web/20171015035606/http://faculty.vassar.edu/lowry/ch8pt1.html

[2]

“卡方检验”,zh.wikipedia.org/wiki/卡方检验

[3]

“G 检验”,[zh.wikipedia.org/wiki/G 检验](https://zh.wikipedia.org/wiki/G 检验)

[4]

Sokal, R. R. 和 Rohlf, F. J. “生物统计学原理与实践”,纽约:Freeman(1981)

[5]

Cressie, N. 和 Read, T. R. C.,“多项式拟合优度检验”,J. Royal Stat. Soc. Series B,Vol. 46, No. 3 (1984),pp. 440-464。

例子

(有关更多示例,请参阅 chisquare。)

当仅提供 f_obs 时,假定期望频率是均匀的,并由观察频率的平均值给出。在这里,我们执行 G 检验(即使用对数似然比统计量):

>>> import numpy as np
>>> from scipy.stats import power_divergence
>>> power_divergence([16, 18, 16, 14, 12, 12], lambda_='log-likelihood')
(2.006573162632538, 0.84823476779463769) 

可以使用 f_exp 参数给出期望频率:

>>> power_divergence([16, 18, 16, 14, 12, 12],
...                  f_exp=[16, 16, 16, 16, 16, 8],
...                  lambda_='log-likelihood')
(3.3281031458963746, 0.6495419288047497) 

f_obs 是二维时,默认情况下,将测试应用于每一列。

>>> obs = np.array([[16, 18, 16, 14, 12, 12], [32, 24, 16, 28, 20, 24]]).T
>>> obs.shape
(6, 2)
>>> power_divergence(obs, lambda_="log-likelihood")
(array([ 2.00657316,  6.77634498]), array([ 0.84823477,  0.23781225])) 

通过设置 axis=None,可以将测试应用于数组中的所有数据,这等效于将测试应用于扁平化的数组。

>>> power_divergence(obs, axis=None)
(23.31034482758621, 0.015975692534127565)
>>> power_divergence(obs.ravel())
(23.31034482758621, 0.015975692534127565) 

ddof 是要对默认自由度进行的更改。

>>> power_divergence([16, 18, 16, 14, 12, 12], ddof=1)
(2.0, 0.73575888234288467) 

通过将测试统计量与 ddof 广播来计算 p 值。

>>> power_divergence([16, 18, 16, 14, 12, 12], ddof=[0,1,2])
(2.0, array([ 0.84914504,  0.73575888,  0.5724067 ])) 

f_obsf_exp 也在广播中使用。在下面的例子中,f_obs 的形状为 (6,),f_exp 的形状为 (2, 6),因此广播 f_obsf_exp 的结果形状为 (2, 6)。要计算所需的卡方统计量,我们必须使用 axis=1

>>> power_divergence([16, 18, 16, 14, 12, 12],
...                  f_exp=[[16, 16, 16, 16, 16, 8],
...                         [8, 20, 20, 16, 12, 12]],
...                  axis=1)
(array([ 3.5 ,  9.25]), array([ 0.62338763,  0.09949846])) 

scipy.stats.ttest_rel

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ttest_rel.html#scipy.stats.ttest_rel

scipy.stats.ttest_rel(a, b, axis=0, nan_policy='propagate', alternative='two-sided', *, keepdims=False)

计算 a 和 b 的两个相关样本的 t 检验。

这是针对两个相关或重复样本具有相同平均(预期)值的零假设的检验。

参数:

a, b类似数组

数组必须具有相同的形状。

axis整数或 None,默认值:0

如果是 int,则是在计算统计量时输入的轴。输入的每个轴切片(例如行)的统计量将出现在输出的相应元素中。如果为None,则在计算统计量之前将展平输入。

nan_policy

定义如何处理输入的 NaN 值。

  • propagate: 如果在计算统计量的轴切片(例如行)中存在 NaN,则输出的相应条目将为 NaN。

  • omit: 在执行计算时将省略 NaN。如果沿着计算统计量的轴切片中剩余的数据不足,输出的相应条目将为 NaN。

  • 如果存在 NaN,则会引发ValueError

alternative,可选

定义备择假设。有以下选项可用(默认为‘two-sided’):

  • ‘two-sided’: 样本基础分布的均值不相等。

  • ‘less’: 第一个样本底层分布的均值小于第二个样本底层分布的均值。

  • ‘greater’: 第一个样本底层分布的均值大于第二个样本底层分布的均值。

从 1.6.0 版本开始。

keepdims布尔型,默认值:False

如果设置为 True,则减少的轴将作为大小为一的维度留在结果中。使用此选项,结果将正确地对输入数组进行广播。

返回:

resultTtestResult

一个带有以下属性的对象:

统计浮点数或数组

t 统计量。

p 值浮点数或数组

与给定备择假设相关的 p 值。

dffloat 或数组

在计算 t 统计量时使用的自由度数量;这比样本的大小少一个(a.shape[axis])。

从 1.10.0 版本开始。

此对象还具有以下方法:

confidence_interval(置信水平=0.95)

为给定置信水平计算群体均值差异的置信区间。置信区间以namedtuple的形式返回,包含lowhigh字段。

从 1.10.0 版本开始。

注意事项

使用示例包括同一组学生在不同考试中的成绩,或者从同一单位重复抽样。该测试评估了平均分数在样本(例如考试)之间是否显著不同。如果观察到一个较大的 p 值,例如大于 0.05 或者 0.1,则我们无法拒绝相同平均分数的零假设。如果 p 值小于阈值,例如 1%、5% 或 10%,则我们拒绝平均值相等的零假设。小的 p 值与大的 t 统计量相关联。

t 统计量计算为 np.mean(a - b)/se,其中 se 是标准误差。因此,当 a - b 的样本均值大于零时,t 统计量为正,当 a - b 的样本均值小于零时,t 统计量为负。

从 SciPy 1.9 开始,np.matrix 输入(不推荐用于新代码)在执行计算前会被转换为 np.ndarray。在这种情况下,输出将是一个适当形状的标量或者 np.ndarray,而不是一个二维的 np.matrix。类似地,虽然掩码数组的掩码元素被忽略,输出将是一个适当形状的标量或者 np.ndarray,而不是具有 mask=False 的掩码数组。

参考文献

en.wikipedia.org/wiki/T-test#Dependent_t-test_for_paired_samples

示例

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng() 
>>> rvs1 = stats.norm.rvs(loc=5, scale=10, size=500, random_state=rng)
>>> rvs2 = (stats.norm.rvs(loc=5, scale=10, size=500, random_state=rng)
...         + stats.norm.rvs(scale=0.2, size=500, random_state=rng))
>>> stats.ttest_rel(rvs1, rvs2)
TtestResult(statistic=-0.4549717054410304, pvalue=0.6493274702088672, df=499)
>>> rvs3 = (stats.norm.rvs(loc=8, scale=10, size=500, random_state=rng)
...         + stats.norm.rvs(scale=0.2, size=500, random_state=rng))
>>> stats.ttest_rel(rvs1, rvs3)
TtestResult(statistic=-5.879467544540889, pvalue=7.540777129099917e-09, df=499) 

scipy.stats.wilcoxon

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.wilcoxon.html#scipy.stats.wilcoxon

scipy.stats.wilcoxon(x, y=None, zero_method='wilcox', correction=False, alternative='two-sided', method='auto', *, axis=0, nan_policy='propagate', keepdims=False)

计算 Wilcoxon 符号秩检验。

Wilcoxon 符号秩检验检验相关配对样本来自相同分布的零假设。特别地,它检验 x - y 的差异分布是否关于零对称。它是配对 T 检验的非参数版本。

参数:

x 类似数组

要么是第一组测量值(在这种情况下 y 是第二组测量值),要么是两组测量值的差异(在这种情况下 y 不应指定)。必须是一维的。

y 类似数组,可选

要么是第二组测量值(如果 x 是第一组测量值),要么未指定(如果 x 是两组测量值之间的差异)。必须是一维的。

警告

当提供 y 时,wilcoxon 根据 d = x - y 的绝对值的排名计算检验统计量。减法中的舍入误差可能导致 d 的元素在确切算术时被分配不同的排名,即使它们会因精确算术而绑定。与分开传递 xy 不同,考虑计算差异 x - y,必要时四舍五入以确保只有真正唯一的元素在数值上是不同的,并将结果作为 x 传递,将 y 保留为默认值(None)。

zero_method {“wilcox”, “pratt”, “zsplit”},可选

处理具有相等值的观测对(“零差异”或“零”的)有不同的约定。

  • “wilcox”:丢弃所有零差异(默认);参见 [4]

  • “pratt”:在排名过程中包括零差异,但删除零的排名(更保守);参见 [3]。在这种情况下,正态近似调整如同 [5]

  • “zsplit”:在排名过程中包括零差异,并将零排名分为正负两部分。

correction 布尔型,可选

如果为 True,在使用正态近似时,通过调整 Wilcoxon 秩统计量向均值调整 0.5 来应用连续性校正。默认为 False。

alternative {“two-sided”, “greater”, “less”},可选

定义备择假设。默认为‘two-sided’。在以下内容中,让 d 表示配对样本之间的差异:如果同时提供 xy,则 d = x - y,否则 d = x

  • ‘two-sided’:d 底层分布不对称于零。

  • ‘less’:d 底层分布在关于零对称的分布上是随机小于的。

  • ‘greater’:d 底层分布在关于零对称的分布上是随机大于的。

method,可选

计算 p 值的方法,请参见备注。默认是“auto”。

axis整数或 None,默认为 0

如果是整数,则沿着输入的轴计算统计量。输入的每个轴切片(例如行)的统计量将出现在输出的相应元素中。如果是 None,则在计算统计量之前会展平输入。

nan_policy

定义如何处理输入的 NaN 值。

  • propagate:如果轴切片(例如行)中存在 NaN,则输出的相应条目将为 NaN。

  • omit:在执行计算时将省略 NaN。如果沿着计算统计量的轴切片中剩余的数据不足,输出的相应条目将为 NaN。

  • raise:如果存在 NaN,则会引发 ValueError

keepdimsbool,默认为 False

如果设置为 True,则被减少的轴将作为大小为一的维度留在结果中。使用此选项,结果将正确地与输入数组进行广播。

返回:

一个具有以下属性的对象。

statistic类似数组

如果 alternative 是“双侧”,则是差异排名之和(无论是正还是负)。否则是正差异的排名之和。

pvalue类似数组

测试的 p 值取决于 alternativemethod

zstatistic类似数组

method = 'approx' 时,这是标准化的 z 统计量:

z = (T - mn - d) / se 

其中 T 是如上定义的 statisticmn 是零假设下分布的均值,d 是连续性校正,se 是标准误差。当 method != 'approx' 时,该属性不可用。

另请参阅

kruskalmannwhitneyu

备注

在以下内容中,让 d 表示成对样本之间的差异:如果提供了 xy,则 d = x - y,否则 d = x。假设所有 d 的元素都是独立同分布的观察值,并且所有元素都是不同且非零的。

  • len(d) 足够大时,标准化检验统计量(如上的 zstatistic)的零分布近似为正态分布,此时可以使用 method = 'approx' 计算 p 值。

  • len(d) 较小时,正态近似可能不准确,推荐使用 method='exact'(尽管执行时间会增加)。

  • 默认情况下,method='auto' 在两者之间选择:当 len(d) <= 50 时,使用精确方法;否则,使用近似方法。

“并列”(即d的所有元素都不唯一)和“零”(即d的元素为零)的存在改变了检验统计量的零分布,当method='exact'时,不再计算精确的 p 值。如果method='approx',则调整了 z 统计量以更准确地与标准正态分布进行比较,但对于有限样本大小,标准正态分布仍然只是 z 统计量真实零分布的近似。关于在零和/或并列存在时,哪种方法最准确地逼近小样本 p 值,参考文献中并无明确共识。无论如何,这是当使用wilcoxonpymethod='auto': ``method='exact'用于len(d) <= 50 并且没有零时的行为;否则,将使用method='approx'

从 SciPy 1.9 开始,不推荐新代码使用的np.matrix输入在计算执行前会转换为np.ndarray。在这种情况下,输出将是合适形状的标量或np.ndarray,而不是二维np.matrix。同样,虽然忽略了掩码数组的掩码元素,输出将是标量或np.ndarray,而不是带有mask=False的掩码数组。

参考文献

[1]

威尔科克森符号秩检验

[2]

Conover, W.J., 实用的非参数统计,1971 年。

[3]

Pratt, J.W., 关于威尔科克森符号秩程序中的零和并列的备注,美国统计协会杂志,第 54 卷,1959 年,第 655-667 页。DOI:10.1080/01621459.1959.10501526

[4] (1,2)

Wilcoxon, F., 通过排名方法进行个体比较,生物统计学通报,第 1 卷,1945 年,第 80-83 页。DOI:10.2307/3001968

[5]

Cureton, E.E., 当零差异存在时,符号秩采样分布的正态近似,美国统计协会杂志,第 62 卷,1967 年,第 1068-1069 页。DOI:10.1080/01621459.1967.10500917

例子

[4]中,自交与异交玉米植株的高度差异如下所示:

>>> d = [6, 8, 14, 16, 23, 24, 28, 29, 41, -48, 49, 56, 60, -67, 75] 

自交植株似乎更高。为了检验没有高度差异的零假设,我们可以应用双侧检验:

>>> from scipy.stats import wilcoxon
>>> res = wilcoxon(d)
>>> res.statistic, res.pvalue
(24.0, 0.041259765625) 

因此,我们会在 5%的置信水平下拒绝零假设,得出两组之间存在高度差异的结论。为了确认差异的中位数可以假定为正,我们使用:

>>> res = wilcoxon(d, alternative='greater')
>>> res.statistic, res.pvalue
(96.0, 0.0206298828125) 

这表明,在 5%的置信水平下,可以拒绝中位数为负的零假设,支持中位数大于零的备择假设。上述 p 值是精确的。使用正态近似得到非常相似的值:

>>> res = wilcoxon(d, method='approx')
>>> res.statistic, res.pvalue
(24.0, 0.04088813291185591) 

在单侧情况下(正差异的秩和),统计量变为 96,而在双侧情况下(零上下秩和的最小值),统计量为 24。

在上述示例中,提供了配对植物高度差异直接给wilcoxon。或者,wilcoxon接受等长的两个样本,计算配对元素之间的差异,然后进行测试。考虑样本 xy

>>> import numpy as np
>>> x = np.array([0.5, 0.825, 0.375, 0.5])
>>> y = np.array([0.525, 0.775, 0.325, 0.55])
>>> res = wilcoxon(x, y, alternative='greater')
>>> res
WilcoxonResult(statistic=5.0, pvalue=0.5625) 

请注意,如果我们手动计算差异,测试结果将会有所不同:

>>> d = [-0.025, 0.05, 0.05, -0.05]
>>> ref = wilcoxon(d, alternative='greater')
>>> ref
WilcoxonResult(statistic=6.0, pvalue=0.4375) 

显著的差异是由于 x-y 结果中的舍入误差造成的:

>>> d - (x-y)
array([2.08166817e-17, 6.93889390e-17, 1.38777878e-17, 4.16333634e-17]) 

即使我们预期所有 (x-y)[1:] 的元素具有相同的幅度 0.05,实际上它们的幅度略有不同,因此在测试中被分配了不同的秩。在执行测试之前,考虑计算 d 并根据需要调整,以确保理论上相同的值在数值上不是不同的。例如:

>>> d2 = np.around(x - y, decimals=3)
>>> wilcoxon(d2, alternative='greater')
WilcoxonResult(statistic=6.0, pvalue=0.4375) 

scipy.stats.linregress

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.linregress.html#scipy.stats.linregress

scipy.stats.linregress(x, y=None, alternative='two-sided')

为两组测量计算线性最小二乘回归。

参数:

x, yarray_like

两组测量值。两个数组应具有相同的长度。如果仅给定 x(并且 y=None),则它必须是一个二维数组,其中一个维度的长度为 2。然后通过沿长度为 2 的维度分割数组来找到两组测量值。在 y=Nonex 是一个 2x2 数组的情况下,linregress(x) 等同于 linregress(x[0], x[1])

alternative, optional

定义备择假设。默认为‘two-sided’。提供以下选项:

  • ‘two-sided’: 回归线的斜率非零

  • ‘less’: 回归线的斜率小于零

  • ‘greater’: 回归线的斜率大于零

自 1.7.0 版本新增。

返回:

resultLinregressResult 实例

返回值是一个带有以下属性的对象:

斜率 float

回归线的斜率。

截距 float

回归线的截距。

rvaluefloat

Pearson 相关系数。rvalue 的平方等于确定系数。

pvaluefloat

用 t-分布的 Wald 检验的检验统计量进行假设检验的 p 值,其零假设是斜率为零。参见上述 alternative 来获取备择假设。

stderrfloat

在残差正态性假设下,估计斜率(梯度)的标准误差。

intercept_stderrfloat

在残差正态性假设下,估计截距的标准误差。

另请参阅

scipy.optimize.curve_fit

使用非线性最小二乘拟合函数到数据。

scipy.optimize.leastsq

最小化一组方程的平方和。

注意事项

将缺失值视为成对处理:如果 x 中的值缺失,则 y 中对应的值被屏蔽。

为了与较早版本的 SciPy 兼容,返回值的行为类似于长度为 5 的 namedtuple,具有字段 slopeinterceptrvaluepvaluestderr,因此可以继续编写:

slope, intercept, r, p, se = linregress(x, y) 

然而,使用该风格时,截距的标准误差不可用。为了访问所有计算值,包括截距的标准误差,请使用返回值作为具有属性的对象,例如:

result = linregress(x, y)
print(result.intercept, result.intercept_stderr) 

示例

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import stats
>>> rng = np.random.default_rng() 

生成一些数据:

>>> x = rng.random(10)
>>> y = 1.6*x + rng.random(10) 

执行线性回归:

>>> res = stats.linregress(x, y) 

决定系数(R 平方):

>>> print(f"R-squared: {res.rvalue**2:.6f}")
R-squared: 0.717533 

将数据与拟合线一起绘制:

>>> plt.plot(x, y, 'o', label='original data')
>>> plt.plot(x, res.intercept + res.slope*x, 'r', label='fitted line')
>>> plt.legend()
>>> plt.show() 

../../_images/scipy-stats-linregress-1_00_00.png

计算斜率和截距的 95% 置信区间:

>>> # Two-sided inverse Students t-distribution
>>> # p - probability, df - degrees of freedom
>>> from scipy.stats import t
>>> tinv = lambda p, df: abs(t.ppf(p/2, df)) 
>>> ts = tinv(0.05, len(x)-2)
>>> print(f"slope (95%): {res.slope:.6f} +/- {ts*res.stderr:.6f}")
slope (95%): 1.453392 +/- 0.743465
>>> print(f"intercept (95%): {res.intercept:.6f}"
...       f" +/- {ts*res.intercept_stderr:.6f}")
intercept (95%): 0.616950 +/- 0.544475 

scipy.stats.pearsonr

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.pearsonr.html#scipy.stats.pearsonr

scipy.stats.pearsonr(x, y, *, alternative='two-sided', method=None)

用于测试非相关性的 Pearson 相关系数和 p 值。

Pearson 相关系数 [1] 用于衡量两个数据集之间的线性关系。与其他相关系数一样,此系数在 -1 到 +1 之间变化,0 表示无相关性。相关系数为 -1 或 +1 表示精确的线性关系。正相关表示随着 x 的增加,y 也增加。负相关表示随着 x 的增加,y 减少。

此函数还执行零假设测试,即样本所代表的分布是不相关和正态分布的。(有关输入非正态对相关系数分布影响的讨论,请参见 Kowalski [3]。)p 值大致指示不相关系统生成具有至少与这些数据集计算的 Pearson 相关性一样极端的数据集的概率。

参数:

x(N,) array_like

输入数组。

y(N,) array_like

输入数组。

alternative,可选

定义备择假设。默认为 ‘two-sided’。可用的选项包括:

  • ‘two-sided’: 相关性非零

  • ‘less’: 相关性为负(小于零)

  • ‘greater’: 相关性为正(大于零)

自版本 1.9.0 新增。

methodResamplingMethod,可选

定义了计算 p 值的方法。如果 methodPermutationMethod/MonteCarloMethod 的实例,则使用提供的配置选项和其他适当的设置使用 scipy.stats.permutation_test/scipy.stats.monte_carlo_test 计算 p 值。否则,按照说明文档计算 p 值。

自版本 1.11.0 新增。

返回:

resultPearsonRResult

一个具有以下属性的对象:

statisticfloat

Pearson 乘积矩相关系数。

pvaluefloat

与选择的备择假设相关的 p 值。

该对象具有以下方法:

confidence_interval(confidence_level, method)

此函数计算给定置信水平的相关系数统计量的置信区间。置信区间以namedtuple的形式返回,字段为lowhigh。如果未提供method,则使用 Fisher 变换计算置信区间[1]。如果methodBootstrapMethod的一个实例,则使用scipy.stats.bootstrap根据提供的配置选项和其他适当的设置计算置信区间。在某些情况下,由于重采样退化,置信限可能为 NaN,这在非常小的样本(~6 个观测值)中很典型。

警告:

ConstantInputWarning

如果输入为常量数组,则引发警告。在这种情况下,相关系数未定义,因此返回np.nan

NearConstantInputWarning

如果输入“几乎”是常量,则引发警告。如果x数组被认为几乎是常量,则norm(x - mean(x)) < 1e-13 * abs(mean(x))。在这种情况下,计算中x - mean(x)的数值误差可能导致 r 的不准确计算。

另见:

spearmanr

Spearman 秩相关系数。

kendalltau

Kendall's tau,用于有序数据的相关度量。

注:

相关系数计算方法如下:

[r = \frac{\sum (x - m_x) (y - m_y)} {\sqrt{\sum (x - m_x)² \sum (y - m_y)²}}]

其中(m_x)为向量 x 的均值,(m_y)为向量 y 的均值。

在假设 x 和 y 来自独立正态分布(因此总体相关系数为 0)的条件下,样本相关系数 r 的概率密度函数为([1], [2]):

[f(r) = \frac{{(1-r²)}^{n/2-2}}{\mathrm{B}(\frac{1}{2},\frac{n}{2}-1)}]

其中 n 为样本数,B 为 Beta 函数。有时称为 r 的确切分布。当method参数保持默认值(None)时,这是用于计算 p 值的分布,例如pearsonr。r 的分布是在区间[-1, 1]上的 Beta 分布,具有相等的形状参数 a = b = n/2 - 1。在 SciPy 的 Beta 分布实现中,r 的分布如下:

dist = scipy.stats.beta(n/2 - 1, n/2 - 1, loc=-1, scale=2) 

pearsonr返回的默认 p 值是双侧 p 值。对于给定相关系数 r 的样本,p 值是随机抽样 x'和 y'来自零相关总体时 abs(r')大于或等于 abs(r)的概率。根据上述dist对象,给定 r 和长度 n 的 p 值可以计算为:

p = 2*dist.cdf(-abs(r)) 

当 n 为 2 时,上述连续分布未定义。可以将 beta 分布在形状参数 a 和 b 接近 a = b = 0 时解释为具有 r = 1 和 r = -1 的离散分布。更直接地,可以观察到,鉴于数据 x = [x1, x2]和 y = [y1, y2],假设 x1 != x2 和 y1 != y2,则 r 的唯一可能值为 1 和-1。因为对于长度为 2 的任意样本 x'和 y',abs(r')始终为 1,所以长度为 2 的样本的双侧 p 值始终为 1。

为了向后兼容,返回的对象也像长度为二的元组,其中保存统计量和 p 值。

参考文献

[1] (1,2,3)

“皮尔逊相关系数”,维基百科,en.wikipedia.org/wiki/Pearson_correlation_coefficient

[2]

Student,“相关系数的可能误差”,生物统计学,第 6 卷,第 2-3 期,1908 年 9 月 1 日,第 302-310 页。

[3]

C. J. Kowalski,“关于非正态对样本积矩相关系数分布的影响” 皇家统计学会杂志。C 系列(应用统计学),第 21 卷,第 1 期(1972 年),第 1-12 页。

例如

>>> import numpy as np
>>> from scipy import stats
>>> x, y = [1, 2, 3, 4, 5, 6, 7], [10, 9, 2.5, 6, 4, 3, 2]
>>> res = stats.pearsonr(x, y)
>>> res
PearsonRResult(statistic=-0.828503883588428, pvalue=0.021280260007523286) 

执行测试的精确排列版本:

>>> rng = np.random.default_rng()
>>> method = stats.PermutationMethod(n_resamples=np.inf, random_state=rng)
>>> stats.pearsonr(x, y, method=method)
PearsonRResult(statistic=-0.828503883588428, pvalue=0.028174603174603175) 

在空假设下执行检验,即数据来自均匀分布:

>>> method = stats.MonteCarloMethod(rvs=(rng.uniform, rng.uniform))
>>> stats.pearsonr(x, y, method=method)
PearsonRResult(statistic=-0.828503883588428, pvalue=0.0188) 

生成渐近 90%置信区间:

>>> res.confidence_interval(confidence_level=0.9)
ConfidenceInterval(low=-0.9644331982722841, high=-0.3460237473272273) 

而对于自举置信区间:

>>> method = stats.BootstrapMethod(method='BCa', random_state=rng)
>>> res.confidence_interval(confidence_level=0.9, method=method)
ConfidenceInterval(low=-0.9983163756488651, high=-0.22771001702132443)  # may vary 

如果 y = a + b*x + e,其中 a,b 是常数,e 是随机误差项,假设与 x 独立。为简单起见,假设 x 是标准正态分布,a=0,b=1,让 e 遵循均值为零,标准差为 s>0 的正态分布。

>>> rng = np.random.default_rng()
>>> s = 0.5
>>> x = stats.norm.rvs(size=500, random_state=rng)
>>> e = stats.norm.rvs(scale=s, size=500, random_state=rng)
>>> y = x + e
>>> stats.pearsonr(x, y).statistic
0.9001942438244763 

这应该接近所给出的确切值。

>>> 1/np.sqrt(1 + s**2)
0.8944271909999159 

对于 s=0.5,我们观察到高度相关性。通常,噪声的大方差会降低相关性,而误差方差接近零时,相关性接近于 1。

需要记住,没有相关性并不意味着独立,除非(x,y)是联合正态的。即使在存在非常简单的依赖结构时,相关性也可能为零:如果 X 服从标准正态分布,则令 y = abs(x)。注意,x 和 y 之间的相关性为零。确实,由于 x 的期望为零,cov(x,y) = E[xy]。根据定义,这等于 E[xabs(x)],由于对称性,这是零。以下代码行说明了这一观察:

>>> y = np.abs(x)
>>> stats.pearsonr(x, y)
PearsonRResult(statistic=-0.05444919272687482, pvalue=0.22422294836207743) 

非零相关系数可能具有误导性。例如,如果 X 符合标准正态分布,定义 y = x 如果 x < 0,否则 y = 0。简单的计算显示 corr(x, y) = sqrt(2/Pi) = 0.797…,表明高度相关:

>>> y = np.where(x < 0, x, 0)
>>> stats.pearsonr(x, y)
PearsonRResult(statistic=0.861985781588, pvalue=4.813432002751103e-149) 

这是不直观的,因为如果我们对 x 和 y 进行抽样,当 x 大于零时,x 和 y 之间没有依赖关系,在大约一半的情况下会发生这种情况。

scipy.stats.spearmanr

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.spearmanr.html#scipy.stats.spearmanr

scipy.stats.spearmanr(a, b=None, axis=0, nan_policy='propagate', alternative='two-sided')

计算 Spearman 相关系数及其相关的 p 值。

Spearman 秩相关系数是两个数据集之间单调关系的非参数测量。与其他相关系数一样,其值在 -1 到 +1 之间,其中 0 表示没有相关性。相关系数为 -1 或 +1 表示精确的单调关系。正相关表示随着 x 的增加,y 也增加。负相关表示随着 x 的增加,y 减少。

p 值大致表示无相关系统生成具有与这些数据集计算出的 Spearman 相关性至少一样极端的数据集的概率。虽然计算 p 值不对样本下面的分布做出强烈的假设,但仅适用于非常大的样本(>500 观测)。对于较小的样本大小,请考虑置换检验(参见下面的示例部分)。

参数:

a, b1D 或 2D array_like,b 是可选的

包含多个变量和观测值的一个或两个 1-D 或 2-D 数组。当这些为 1-D 时,每个表示单个变量的观测值向量。在 2-D 情况下的行为,请参见下面的axis。两个数组在axis维度上需要具有相同的长度。

axisint 或 None,可选

如果axis=0(默认),则每列代表一个变量,行中包含观测值。如果axis=1,则关系被转置:每行表示一个变量,而列包含观测值。如果axis=None,则两个数组都会被展平。

nan_policy,可选

定义了如何处理输入包含 NaN 的情况。以下选项可用(默认为‘propagate’):

  • ‘propagate’:返回 NaN

  • ‘raise’:抛出错误

  • ‘omit’:执行计算时忽略 NaN 值

alternative,可选

定义了备择假设。默认为‘two-sided’。以下选项可用:

  • ‘two-sided’:相关性为非零

  • ‘less’:相关性为负(小于零)

  • ‘greater’:相关性为正(大于零)

新版本为 1.7.0。

返回:

resSignificanceResult

一个包含属性的对象:

statisticfloat 或 ndarray(2-D 方阵)

Spearman 相关系数矩阵或相关系数(如果仅给出 2 个变量作为参数)。相关系数矩阵是方阵,其长度等于ab组合后的总变量数(列或行)。

pvaluefloat

p 值用于一个假设检验,其零假设是两个样本没有顺序相关性。参见上面的alternative用于备择假设。pvalue具有与statistic相同的形状。

警告:

ConstantInputWarning

如果输入是一个常数数组,则引发此警告。在这种情况下,相关系数未定义,因此返回np.nan

参考文献

[1]

Zwillinger, D. 和 Kokoska, S. (2000). CRC 标准概率和统计表格与公式. Chapman & Hall: New York. 2000. Section 14.7

[2]

Kendall, M. G. 和 Stuart, A. (1973). 统计学的高级理论,卷 2:推理与关系. Griffin. 1973. Section 31.18

[3]

Kershenobich, D., Fierro, F. J., & Rojkind, M. (1970). 游离脯氨酸与人类肝硬化中胶原含量的关系. The Journal of Clinical Investigation, 49(12), 2246-2249.

[4]

Hollander, M., Wolfe, D. A., & Chicken, E. (2013). 非参数统计方法. John Wiley & Sons.

[5]

B. Phipson 和 G. K. Smyth. “置换 P 值永远不应为零:当置换随机抽取时计算精确 P 值。” 遗传和分子生物统计应用 9.1 (2010).

[6]

Ludbrook, J., & Dudley, H. (1998). 为什么在生物医学研究中置换测试优于 t 和 F 测试. The American Statistician, 52(2), 127-132.

示例

考虑以下来自[3]的数据,研究了不健康人类肝脏中游离脯氨酸(一种氨基酸)和总胶原(经常存在于结缔组织中的蛋白质)之间的关系。

下面的xy数组记录了这两种化合物的测量值。这些观察值是成对的:每个游离脯氨酸测量是在相同的肝脏中以相同的索引进行的总胶原测量。

>>> import numpy as np
>>> # total collagen (mg/g dry weight of liver)
>>> x = np.array([7.1, 7.1, 7.2, 8.3, 9.4, 10.5, 11.4])
>>> # free proline (μ mole/g dry weight of liver)
>>> y = np.array([2.8, 2.9, 2.8, 2.6, 3.5, 4.6, 5.0]) 

这些数据在[4]中使用斯皮尔曼相关系数进行了分析,这是一种对样本间单调相关性敏感的统计量。

>>> from scipy import stats
>>> res = stats.spearmanr(x, y)
>>> res.statistic
0.7000000000000001 

这一统计量的值在样本间具有强烈正序相关性时趋向于高(接近 1),在样本间具有强烈负序相关性时趋向于低(接近-1),对于弱序相关性的样本,其大小接近于零。

该测试通过将统计量的观察值与空假设的空分布进行比较来进行。在空假设下,总胶原和游离脯氨酸测量是独立的。

对于这个测试,统计量可以转换,使得大样本的空假设分布为自由度为len(x) - 2的学生 t 分布。

>>> import matplotlib.pyplot as plt
>>> dof = len(x)-2  # len(x) == len(y)
>>> dist = stats.t(df=dof)
>>> t_vals = np.linspace(-5, 5, 100)
>>> pdf = dist.pdf(t_vals)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def plot(ax):  # we'll reuse this
...     ax.plot(t_vals, pdf)
...     ax.set_title("Spearman's Rho Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> plot(ax)
>>> plt.show() 

../../_images/scipy-stats-spearmanr-1_00_00.png

比较通过 p 值来量化:在两侧检验中,统计量为正数时,零分布中大于变换统计量的元素和零分布中小于观测统计量的负值都被视为“更极端”。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> rs = res.statistic  # original statistic
>>> transformed = rs * np.sqrt(dof / ((rs+1.0)*(1.0-rs)))
>>> pvalue = dist.cdf(-transformed) + dist.sf(transformed)
>>> annotation = (f'p-value={pvalue:.4f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (2.7, 0.025), (3, 0.03), arrowprops=props)
>>> i = t_vals >= transformed
>>> ax.fill_between(t_vals[i], y1=0, y2=pdf[i], color='C0')
>>> i = t_vals <= -transformed
>>> ax.fill_between(t_vals[i], y1=0, y2=pdf[i], color='C0')
>>> ax.set_xlim(-5, 5)
>>> ax.set_ylim(0, 0.1)
>>> plt.show() 

../../_images/scipy-stats-spearmanr-1_01_00.png

>>> res.pvalue
0.07991669030889909  # two-sided p-value 

如果 p 值“小” - 也就是说,从独立分布中抽样产生这样极端统计量值的概率很低 - 这可能被视为反对零假设,赞同替代假设:总胶原蛋白和游离脯氨酸的分布独立。请注意:

  • 反之则不成立;也就是说,该检验不用于提供零假设的证据。

  • 被视为“小”的值的阈值是在分析数据之前作出的选择,考虑到假阳性(错误拒绝零假设)和假阴性(未能拒绝假零假设)的风险[5]

  • 小的 p 值不是大效应的证据;而是只能提供“显著”效应的证据,意味着在零假设下发生这样极端值的概率很低。

假设在执行实验之前,作者有理由预测总胶原蛋白和游离脯氨酸测量之间存在正相关,并选择评估零假设对单侧替代的合理性:游离脯氨酸与总胶原蛋白呈正序相关。在这种情况下,只有零分布中那些与观察统计量一样大或更大的值被认为更加极端。

>>> res = stats.spearmanr(x, y, alternative='greater')
>>> res.statistic
0.7000000000000001  # same statistic
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> pvalue = dist.sf(transformed)
>>> annotation = (f'p-value={pvalue:.6f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (3, 0.018), (3.5, 0.03), arrowprops=props)
>>> i = t_vals >= transformed
>>> ax.fill_between(t_vals[i], y1=0, y2=pdf[i], color='C0')
>>> ax.set_xlim(1, 5)
>>> ax.set_ylim(0, 0.1)
>>> plt.show() 

../../_images/scipy-stats-spearmanr-1_02_00.png

>>> res.pvalue
0.03995834515444954  # one-sided p-value; half of the two-sided p-value 

注意,t 分布提供了零分布的渐近近似;仅对观测值多的样本准确。对于小样本,执行置换检验可能更合适:在总胶原蛋白和游离脯氨酸独立的零假设下,每个游离脯氨酸测量可能与任何总胶原蛋白测量一起被观测到。因此,我们可以通过计算在xy之间每一对元素的统计量来形成一个精确的零分布。

>>> def statistic(x):  # explore all possible pairings by permuting `x`
...     rs = stats.spearmanr(x, y).statistic  # ignore pvalue
...     transformed = rs * np.sqrt(dof / ((rs+1.0)*(1.0-rs)))
...     return transformed
>>> ref = stats.permutation_test((x,), statistic, alternative='greater',
...                              permutation_type='pairings')
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> ax.hist(ref.null_distribution, np.linspace(-5, 5, 26),
...         density=True)
>>> ax.legend(['aymptotic approximation\n(many observations)',
...            f'exact \n({len(ref.null_distribution)} permutations)'])
>>> plt.show() 

../../_images/scipy-stats-spearmanr-1_03_00.png

>>> ref.pvalue
0.04563492063492063  # exact one-sided p-value 

scipy.stats.pointbiserialr

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.pointbiserialr.html#scipy.stats.pointbiserialr

scipy.stats.pointbiserialr(x, y)

计算点二列相关系数及其 p 值。

点二列相关用于衡量二进制变量 x 与连续变量 y 之间的关系。与其他相关系数一样,其取值介于 -1 到 +1 之间,0 表示无相关。相关系数为 -1 或 +1 表示决定性关系。

可以使用快捷公式计算此函数,但结果与 pearsonr 相同。

参数:

xarray_like of bools

输入数组。

yarray_like

输入数组。

返回值:

res: SignificanceResult

一个包含以下属性的对象:

statisticfloat

R 值。

pvaluefloat

双侧 p 值。

注意事项

pointbiserialr 使用具有 n-1 自由度的 t 检验。它相当于 pearsonr

点二列相关的值可以从以下公式计算得出:

[r_{pb} = \frac{\overline{Y_1} - \overline{Y_0}} {s_y} \sqrt{\frac{N_0 N_1} {N (N - 1)}}]

其中 (\overline{Y_{0}}) 和 (\overline{Y_{1}}) 分别是编码为 0 和 1 的度量观测值的均值;(N_{0}) 和 (N_{1}) 分别是编码为 0 和 1 的观测数量;(N) 是所有观测值的总数,(s_{y}) 是所有度量观测值的标准差。

当 (r_{pb}) 的值显著不为零时,完全等同于两组之间均值的显著差异。因此,可以使用具有 (N-2) 自由度的独立组 t 检验来检验 (r_{pb}) 是否为非零。比较两个独立组的 t 统计量与 (r_{pb}) 之间的关系如下:

[t = \sqrt{N - 2}\frac{r_{pb}}{\sqrt{1 - r^{2}_{pb}}}]

参考文献

[1]

J. Lev,“点二列相关系数”,Ann. Math. Statist.,Vol. 20,no.1,pp. 125-126,1949 年。

[2]

R.F. Tate,“离散和连续变量之间的相关性。点二列相关。”,Ann. Math. Statist.,Vol. 25,np. 3,pp. 603-607,1954 年。

[3]

D. Kornbrot,“点二列相关”,载于 Wiley StatsRef:统计参考在线版(eds N. Balakrishnan 等),2014 年。DOI:10.1002/9781118445112.stat06227

例子

>>> import numpy as np
>>> from scipy import stats
>>> a = np.array([0, 0, 0, 1, 1, 1, 1])
>>> b = np.arange(7)
>>> stats.pointbiserialr(a, b)
(0.8660254037844386, 0.011724811003954652)
>>> stats.pearsonr(a, b)
(0.86602540378443871, 0.011724811003954626)
>>> np.corrcoef(a, b)
array([[ 1\.       ,  0.8660254],
 [ 0.8660254,  1\.       ]]) 

scipy.stats.kendalltau

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.kendalltau.html#scipy.stats.kendalltau

scipy.stats.kendalltau(x, y, *, initial_lexsort=<object object>, nan_policy='propagate', method='auto', variant='b', alternative='two-sided')

计算 Kendall’s tau,用于序数数据的相关性测量。

Kendall’s tau 是两个排名之间一致性的度量。接近 1 的值表示强烈一致,接近-1 的值表示强烈不一致。此实现了 Kendall’s tau 的两个变体:tau-b(默认)和 tau-c(也称为 Stuart’s tau-c)。它们唯一的区别在于它们如何被归一化到-1 到 1 的范围内;假设检验(它们的 p 值)是相同的。Kendall’s 最初的 tau-a 没有单独实现,因为在没有并列值的情况下,tau-b 和 tau-c 都归约为 tau-a。

参数:

x, yarray_like

排名数组,形状相同。如果数组不是 1-D,则将其展平为 1-D。

initial_lexsortbool,可选,已弃用

此参数未使用。

自版本 1.10.0 起弃用:kendalltau 关键字参数 initial_lexsort 已弃用,因为未使用且将在 SciPy 1.14.0 中移除。

nan_policy,可选

定义当输入包含 NaN 时如何处理。可选的选项如下(默认为‘propagate’):

  • ‘propagate’:返回 NaN
  • ‘raise’:引发错误
  • ‘omit’:执行计算时忽略 NaN 值

method,可选

定义用于计算 p 值的方法 [5]。可选的选项如下(默认为‘auto’):

  • ‘auto’:根据速度和精度之间的平衡选择适当的方法
  • ‘asymptotic’:对大样本有效的正态近似
  • ‘exact’:计算精确的 p 值,但只能在没有并列值的情况下使用。随着样本量的增加,‘exact’ 计算时间可能会增加,并且结果可能会失去一些精度。

variant,可选

定义返回的 Kendall’s tau 的变体。默认为‘b’。

alternative,可选

定义备择假设。默认为‘two-sided’。可选的选项如下:

  • ‘two-sided’:排名相关性非零

  • ‘less’:排名相关性为负(小于零)

  • ‘greater’: 排名相关性为正(大于零)

返回:

resSignificanceResult

一个对象,包含以下属性:

统计量 float

tau 统计量。

p 值 float

用于假设检验的 p 值,其零假设为无关联,tau = 0。

另请参阅

spearmanr

计算 Spearman 秩相关系数。

theilslopes

计算一组点(x, y)的 Theil-Sen 估计量。

weightedtau

计算 Kendall’s tau 的加权版本。

所使用的 Kendall’s tau 定义是 [2]:

tau_b = (P - Q) / sqrt((P + Q + T) * (P + Q + U))

tau_c = 2 (P - Q) / (n**2 * (m - 1) / m) 

其中 P 是协调对的数量,Q 是不协调对的数量,T 是仅在x中的绑定数,U 是仅在y中的绑定数。如果同一对在xy中都有绑定,则不会添加到 T 或 U 中。n 是样本的总数,m 是在xy中较小的唯一值的数量。

参考文献

[1]

Maurice G. Kendall, “排名相关性的新测量”, Biometrika Vol. 30, No. 1/2, pp. 81-93, 1938.

[2]

Maurice G. Kendall, “在排名问题中处理绑定的方法”, Biometrika Vol. 33, No. 3, pp. 239-251. 1945.

[3]

Gottfried E. Noether, “非参数统计要素”, John Wiley & Sons, 1967.

[4]

Peter M. Fenwick, “用于累积频率表的新数据结构”, Software: Practice and Experience, Vol. 24, No. 3, pp. 327-336, 1994.

[5]

Maurice G. Kendall, “排名相关性方法” (第 4 版), Charles Griffin & Co., 1970.

[6]

Kershenobich, D., Fierro, F. J., & Rojkind, M. (1970). 自由脯氨酸的自由池与人类肝硬化中的胶原含量之间的关系。 The Journal of Clinical Investigation, 49(12), 2246-2249.

[7]

Hollander, M., Wolfe, D. A., & Chicken, E. (2013). 非参数统计方法。 John Wiley & Sons.

[8]

B. Phipson 和 G. K. Smyth. “当置换随机抽取时,置换 P 值永远不应该为零:计算确切的 P 值。” Statistical Applications in Genetics and Molecular Biology 9.1 (2010).

示例

请考虑来自 [6] 的以下数据,该研究了不健康的人类肝脏中自由脯氨酸(一种氨基酸)与总胶原(一种经常在结缔组织中找到的蛋白质)之间的关系。

下面的xy数组记录了两种化合物的测量结果。观察结果是成对的:每个自由脯氨酸的测量值都是从相同的肝脏中取得的,与相同索引处的总胶原测量值对应。

>>> import numpy as np
>>> # total collagen (mg/g dry weight of liver)
>>> x = np.array([7.1, 7.1, 7.2, 8.3, 9.4, 10.5, 11.4])
>>> # free proline (μ mole/g dry weight of liver)
>>> y = np.array([2.8, 2.9, 2.8, 2.6, 3.5, 4.6, 5.0]) 

这些数据在 [7] 中使用了 Spearman’s 相关系数进行分析,这是一种与 Kendall’s tau 类似的统计量,同样对样本之间的序数相关性敏感。让我们使用 Kendall’s tau 进行类似的研究。

>>> from scipy import stats
>>> res = stats.kendalltau(x, y)
>>> res.statistic
0.5499999999999999 

对于具有强正序数相关的样本,该统计量的值往往很高(接近 1),对于具有强负序数相关的样本,该值很低(接近-1),对于具有弱序数相关的样本,该值的数量级很小(接近零)。

通过将统计量的观察值与空假设下的空分布进行比较来执行测试:总胶原和自由脯氨酸测量是独立的空假设的统计量分布。

对于此检验,大样本且无绑定的零分布被近似为具有方差 (2*(2*n + 5))/(9*n*(n - 1)) 的正态分布,其中 n = len(x)

>>> import matplotlib.pyplot as plt
>>> n = len(x)  # len(x) == len(y)
>>> var = (2*(2*n + 5))/(9*n*(n - 1))
>>> dist = stats.norm(scale=np.sqrt(var))
>>> z_vals = np.linspace(-1.25, 1.25, 100)
>>> pdf = dist.pdf(z_vals)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def plot(ax):  # we'll reuse this
...     ax.plot(z_vals, pdf)
...     ax.set_title("Kendall Tau Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> plot(ax)
>>> plt.show() 

../../_images/scipy-stats-kendalltau-1_00_00.png

比较通过 p 值量化:在双侧检验中,统计量为正时,零分布中大于转换后的统计量的值和零分布中小于观察统计量的负值都被认为是“更极端”的。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> pvalue = dist.cdf(-res.statistic) + dist.sf(res.statistic)
>>> annotation = (f'p-value={pvalue:.4f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (0.65, 0.15), (0.8, 0.3), arrowprops=props)
>>> i = z_vals >= res.statistic
>>> ax.fill_between(z_vals[i], y1=0, y2=pdf[i], color='C0')
>>> i = z_vals <= -res.statistic
>>> ax.fill_between(z_vals[i], y1=0, y2=pdf[i], color='C0')
>>> ax.set_xlim(-1.25, 1.25)
>>> ax.set_ylim(0, 0.5)
>>> plt.show() 

../../_images/scipy-stats-kendalltau-1_01_00.png

>>> res.pvalue
0.09108705741631495  # approximate p-value 

请注意,曲线的阴影区域与kendalltau返回的 p 值之间存在轻微差异。这是因为我们的数据存在绑定,并且我们忽略了kendalltau执行的零分布方差的绑定修正。对于没有绑定的样本,我们的图表的阴影区域和kendalltau返回的 p 值会完全匹配。

如果 p 值“小” - 即从独立分布中抽取产生这样一个极端统计量值的概率很低 - 这可能被视为反对零假设的证据,支持备择假设:总胶原蛋白和游离脯氨酸的分布独立。请注意:

  • 反之则不成立;也就是说,该检验不用于提供支持零假设的证据。

  • 被视为“小”的值的阈值是在分析数据之前应该做出的选择,考虑到误报(错误拒绝零假设)和误放大(未能拒绝虚假零假设)的风险[8]

  • 小的 p 值并不支持效应的证据;相反,它们只能提供“显著”效应的证据,即它们在零假设下发生的可能性很低。

对于中等规模无绑定样本,kendalltau 可以精确计算 p 值。然而,在存在绑定的情况下,kendalltau 将采用渐近逼近法。尽管如此,我们可以使用置换检验来精确计算零分布:在总胶原蛋白和游离脯氨酸独立的零假设下,每个游离脯氨酸测量值与任何总胶原蛋白测量值一样可能被观察到。因此,我们可以通过计算在 xy 之间每个可能配对的元素下的统计量来形成一个精确的零分布。

>>> def statistic(x):  # explore all possible pairings by permuting `x`
...     return stats.kendalltau(x, y).statistic  # ignore pvalue
>>> ref = stats.permutation_test((x,), statistic,
...                              permutation_type='pairings')
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> bins = np.linspace(-1.25, 1.25, 25)
>>> ax.hist(ref.null_distribution, bins=bins, density=True)
>>> ax.legend(['aymptotic approximation\n(many observations)',
...            'exact null distribution'])
>>> plot(ax)
>>> plt.show() 

../../_images/scipy-stats-kendalltau-1_02_00.png

>>> ref.pvalue
0.12222222222222222  # exact p-value 

注意,这里计算得到的精确 p 值与上述kendalltau返回的近似值存在显著差异。对于具有绑定的小样本,请考虑执行置换检验以获得更准确的结果。

scipy.stats.weightedtau

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.weightedtau.html#scipy.stats.weightedtau

scipy.stats.weightedtau(x, y, rank=True, weigher=None, additive=True)

计算 Kendall 的加权版本 (\tau)。

加权 (\tau) 是 Kendall (\tau) 的加权版本,在此版本中,高权重的交换比低权重的交换更具影响力。默认参数计算指数加法版本的指数,(\tau_\mathrm h),已被证明在重要和不重要元素之间提供了最佳平衡[1]

加权是通过一个等级数组和一个称重函数定义的,该函数为每个元素分配基于等级的权重(较重要的等级与较小的值相关联,例如,0 是最高可能的等级),然后交换的权重是交换元素等级的权重的和或乘积。默认参数计算 (\tau_\mathrm h):在等级为 (r) 和 (s)(从零开始)的元素之间的交换的权重为 (1/(r+1) + 1/(s+1))。

指定等级数组只有在您有一个外部重要性标准的情况下才有意义。如果像通常发生的那样,您没有一个具体的等级标准在脑海中,那么加权 (\tau) 就是通过使用 (x, y) 和 (y, x) 递减字典序排名得到的值的平均值来定义的。这是默认参数的行为。请注意,这里用于排名的约定(较低的值意味着更高的重要性)与其他 SciPy 统计函数使用的约定相反。

参数:

x, y数组样本

得分数组,形状相同。如果数组不是 1-D,则将其展平为 1-D。

rank整数数组或布尔值的数组,可选

给每个元素分配一个非负的等级。如果为 None,则将使用递减字典序排名 (x, y):更高等级的元素将是具有更大 x 值的元素,使用 y 值来打破并列(特别地,交换 xy 将产生不同的结果)。如果为 False,则将直接使用元素索引作为等级。默认为 True,此时该函数返回使用 (x, y) 和 (y, x) 递减字典序排名得到的值的平均值。

weigher可调用对象,可选

该称重函数必须将非负整数(零表示最重要的元素)映射到非负权重。默认情况下,None 提供双曲线加权,即,排名 (r) 被映射到权重 (1/(r+1))。

additive布尔值,可选

如果为 True,则交换的权重通过添加交换元素的等级的权重来计算;否则,权重将相乘。默认为 True。

返回:

res: SignificanceResult

包含属性的对象:

statisticfloat

加权的τ相关指数。

pvaluefloat

目前为np.nan,因为统计量的空分布未知(即使在加性双曲线情况下也是如此)。

另见

kendalltau

计算 Kendall's tau。

spearmanr

计算 Spearman 等级相关系数。

theilslopes

计算一组点(x,y)的 Theil-Sen 估计器。

注意

此函数使用基于(O(n \log n))的归并排序算法[1],这是肯德尔τ的 Knight 算法的加权扩展[2]。它可以通过将additiverank设置为 False 来计算 Shieh 的加权τ[3],用于排名之间无并列(即排列)的情况,因为[1]中给出的定义是 Shieh 的一般化。

NaNs 被认为是最小可能的分数。

0.19.0 版中的新功能。

参考文献

[1] (1,2,3)

Sebastiano Vigna,《带有并列的排名的加权相关指数》,《第 24 届国际万维网会议论文集》,第 1166-1176 页,ACM,2015 年。

[2]

W.R. Knight,《一种计算 Kendall's Tau 的计算机方法,适用于非分组数据》,《美国统计协会杂志》,第 61 卷,第 314 号,第一部分,第 436-439 页,1966 年。

[3]

Grace S. Shieh,《加权的肯德尔τ统计量》,《统计与概率信函》,第 39 卷,第 1 期,第 17-24 页,1998 年。

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = [12, 2, 1, 12, 2]
>>> y = [1, 4, 7, 1, 0]
>>> res = stats.weightedtau(x, y)
>>> res.statistic
-0.56694968153682723
>>> res.pvalue
nan
>>> res = stats.weightedtau(x, y, additive=False)
>>> res.statistic
-0.62205716951801038 

NaNs 被认为是最小可能的分数:

>>> x = [12, 2, 1, 12, 2]
>>> y = [1, 4, 7, 1, np.nan]
>>> res = stats.weightedtau(x, y)
>>> res.statistic
-0.56694968153682723 

这恰好是 Kendall's tau:

>>> x = [12, 2, 1, 12, 2]
>>> y = [1, 4, 7, 1, 0]
>>> res = stats.weightedtau(x, y, weigher=lambda x: 1)
>>> res.statistic
-0.47140452079103173 
>>> x = [12, 2, 1, 12, 2]
>>> y = [1, 4, 7, 1, 0]
>>> stats.weightedtau(x, y, rank=None)
SignificanceResult(statistic=-0.4157652301037516, pvalue=nan)
>>> stats.weightedtau(y, x, rank=None)
SignificanceResult(statistic=-0.7181341329699028, pvalue=nan) 

scipy.stats.somersd

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.somersd.html#scipy.stats.somersd

scipy.stats.somersd(x, y=None, alternative='two-sided')

计算 Somers' D,一种有序关联的非对称度量。

像 Kendall's (\tau) 一样,Somers' (D) 是两个排名之间对应的一种度量。这两个统计量都考虑了两个排名 (X) 和 (Y) 中协调和不协调对的差异,并且都被归一化,使得接近 1 的值表示强烈一致,接近-1 的值表示强烈不一致。它们在归一化方式上有所不同。为了显示关系,Somers' (D) 可以用 Kendall's (\tau_a) 定义:

[D(Y|X) = \frac{\tau_a(X, Y)}{\tau_a(X, X)}]

假设第一个排名 (X) 有 (r) 个不同的排名,第二个排名 (Y) 有 (s) 个不同的排名。这两个由 (n) 个排名组成的列表也可以看作是一个 (r \times s) 的列联表,其中元素 (i, j) 是排名 (X) 中排名 (i) 和排名 (Y) 中排名 (j) 的对数。因此,somersd 还允许将输入数据提供为单个的二维列联表,而不是两个分开的一维排名。

注意,Somers' (D) 的定义是非对称的:一般来说,(D(Y|X) \neq D(X|Y))。somersd(x, y) 计算的是 Somers' (D(Y|X)):将“行”变量 (X) 视为独立变量,“列”变量 (Y) 视为依赖变量。要计算 Somers' (D(X|Y)),请交换输入列表或转置输入表。

参数:

xarray_like

1D 排名数组,被视为(行)独立变量。或者,一个二维列联表。

yarray_like, optional

如果 x 是一个一维排名数组,y 是相同长度的一维排名数组,被视为(列)依赖变量。如果 x 是二维的,则忽略 y

alternative, optional

定义备择假设。默认为 'two-sided'。可用的选项包括:* 'two-sided': 排名相关不为零 * 'less': 排名相关为负(小于零) * 'greater': 排名相关为正(大于零)

返回:

resSomersDResult

一个 SomersDResult 对象,具有以下字段:

statisticfloat

Somers' (D) 统计量。

pvaluefloat

假设检验的 p 值,其零假设是没有关联,即 (D=0)。更多信息请参阅注释。

table2D array

由排名 xy 形成的列联表(或者如果 x 是二维数组,则为提供的列联表)

参见

kendalltau

计算 Kendall's tau,另一种相关度量。

weightedtau

计算 Kendall’s tau 的加权版本。

spearmanr

计算 Spearman 等级相关系数。

pearsonr

计算 Pearson 相关系数。

注意事项

此函数遵循 [2][3] 的列联表方法。p-值是基于在零假设 (D=0) 下的检验统计分布的渐近逼近计算的。

理论上,基于 Kendall’s (tau) 和 Somers’ (D) 的假设检验应该是相同的。然而,kendalltau 返回的 p-值基于 (X) 和 (Y) 之间独立性的零假设(即从中抽取 (X) 和 (Y) 对的总体包含所有可能对的等数量),这比此处使用的 (D=0) 的零假设更为具体。如果需要独立性的零假设,则可以使用 kendalltau 返回的 p-值和 somersd 返回的统计量,反之亦然。更多信息,请参阅 [2]

按照 SAS 和 R 使用的约定格式化列联表:第一个提供的排名(x)是“行”变量,第二个提供的排名(y)是“列”变量。这与 Somers 的原始论文的约定相反 [1]

参考文献

[1]

Robert H. Somers,《用于序数变量的新的非对称关联度量》,《美国社会学评论》,第 27 卷,第 6 期,799–811 页,1962 年。

[2] (1,2)

Morton B. Brown 和 Jacqueline K. Benedetti,《在二维列联表中检验相关性的抽样行为》,《美国统计协会期刊》第 72 卷,第 358 期,309–315 页,1977 年。

[3]

SAS Institute, Inc.,《频数程序(书摘)》,《SAS/STAT 9.2 用户指南,第二版》,SAS Publishing,2009 年。

[4]

Laerd 统计,《使用 SPSS 统计的 Somers’ d》,《SPSS 统计教程和统计指南》,statistics.laerd.com/spss-tutorials/somers-d-using-spss-statistics.php,访问日期为 2020 年 7 月 31 日。

例子

我们为[4]中的示例计算 Somers' D,其中一位酒店连锁店主想要确定酒店房间清洁度与客户满意度之间的关联。自变量酒店房间清洁度在一个有序尺度上进行排名:“低于平均(1)”,“平均(2)”或“高于平均(3)”。因变量客户满意度在第二个尺度上进行排名:“非常不满意(1)”,“中度不满意(2)”,“既不满意也不满意(3)”,“中度满意(4)”,或“非常满意(5)”。共有 189 位顾客参与了调查,结果转化为一个以酒店房间清洁度为“行”变量和客户满意度为“列”变量的列联表。

27 25 14 7 0
7 14 18 35 12
1 3 2 7 17

例如,27 位顾客将其房间的清洁度排名为“低于平均(1)”,相应的满意度为“非常不满意(1)”。我们按以下方式进行分析。

>>> from scipy.stats import somersd
>>> table = [[27, 25, 14, 7, 0], [7, 14, 18, 35, 12], [1, 3, 2, 7, 17]]
>>> res = somersd(table)
>>> res.statistic
0.6032766111513396
>>> res.pvalue
1.0007091191074533e-27 

Somers' D 统计量的值约为 0.6,表明样本中房间清洁度与客户满意度之间存在正相关关系。 p-value 非常小,表明在零假设下观察到该统计量极端值的概率非常小(我们的样本来自 189 位顾客,整体人群的统计量为零假设)。这支持备择假设,即人群的真实 Somers' D 值不为零。

scipy.stats.siegelslopes

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.siegelslopes.html#scipy.stats.siegelslopes

scipy.stats.siegelslopes(y, x=None, method='hierarchical')

对于点集合(x, y),计算Siegel 估计量

siegelslopes 实现了使用重复中位数进行鲁棒线性回归的方法(参见[1]),以拟合点集(x, y)的直线。该方法对异常值具有 50%的渐近破坏点。

参数:

y 数组型

因变量。

x 数组型或 None,可选

自变量。如果为 None,则使用arange(len(y))代替。

方法

如果是‘层次化’,使用估计的斜率slope估计截距(默认选项)。如果是‘分离’,独立估计截距。详见注释。

返回:

result SiegelslopesResult 实例

返回值是一个具有以下属性的对象:

斜率浮点数

回归线斜率的估计。

截距浮点数

回归线截距的估计。

另请参阅

theilslopes

一种类似的技术,但没有重复中位数

注释

对于n = len(y),将m_j计算为从点(x[j], y[j])到所有其他n-1点的斜率的中位数。然后slope是所有斜率m_j的中位数。可以通过参数method选择两种估计截距的方法。层次化方法使用估计的斜率slope,计算intercept作为y - slope*x的中位数。另一种方法独立估计截距如下:对于每个点(x[j], y[j]),计算通过其余点的所有n-1条线的截距i_j的中位数。intercepti_j的中位数。

该实现计算大小为n的向量的中位数n次,对于大向量可能较慢。有更高效的算法(参见[2]),此处未实现。

为了与 SciPy 旧版本兼容,返回值行为类似于长度为 2 的namedtuple,包含字段slopeintercept,因此可以继续写:

slope, intercept = siegelslopes(y, x) 

参考文献

[1] (1,2)

A. Siegel,“使用重复中位数的鲁棒回归”,Biometrika,Vol. 69,pp. 242-244,1982 年。

[2]

A. Stein 和 M. Werman,“寻找重复中位数回归线”,第三届 ACM-SIAM 离散算法年会论文集,pp. 409-413,1992 年。

示例

>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt 
>>> x = np.linspace(-5, 5, num=150)
>>> y = x + np.random.normal(size=x.size)
>>> y[11:15] += 10  # add outliers
>>> y[-5:] -= 7 

计算斜率和截距。为了比较,还可以使用linregress计算最小二乘拟合:

>>> res = stats.siegelslopes(y, x)
>>> lsq_res = stats.linregress(x, y) 

绘制结果。Siegel 回归线以红色显示。绿色线显示最小二乘拟合以供比较。

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(x, y, 'b.')
>>> ax.plot(x, res[1] + res[0] * x, 'r-')
>>> ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'g-')
>>> plt.show() 

../../_images/scipy-stats-siegelslopes-1.png

scipy.stats.theilslopes

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.theilslopes.html#scipy.stats.theilslopes

scipy.stats.theilslopes(y, x=None, alpha=0.95, method='separate')

计算一组点(x, y)的 Theil-Sen 估计器。

theilslopes 实现了一种鲁棒线性回归的方法。它计算斜率作为所有配对值之间斜率的中位数。

参数:

yarray_like

因变量。

xarray_like 或 None,可选

自变量。如果为 None,则使用arange(len(y))

alphafloat,可选

置信度在 0 到 1 之间,默认为 95% 置信度。请注意,alpha 对称地围绕 0.5,即 0.1 和 0.9 都被解释为“查找 90% 置信区间”。

方法,可选

用于计算截距估计的方法。支持以下方法,

  • ‘joint’: 使用 np.median(y - slope * x) 作为截距。
  • ‘separate’: 使用 np.median(y) - slope * np.median(x)
  • 作为截距。

默认值为‘separate’。

版本 1.8.0 中的新功能。

返回:

resultTheilslopesResult 实例

返回值是一个具有以下属性的对象:

斜率 float

Theil 斜率。

截距 float

Theil 线的截距。

低斜率 float

斜率置信区间的下限

高斜率 float

斜率置信区间的上限

参见

siegelslopes

使用重复中位数的类似技术

注意事项

theilslopes 的实现遵循 [1]。在 [1] 中未定义截距,在这里定义为 median(y) - slope*median(x),这在 [3] 中给出。文献中也有其他截距的定义,例如在 [4] 中的 median(y - slope*x)。确定如何计算截距可以通过参数 method 来确定。由于文献中未涉及,因此没有给出截距的置信区间。

为了与 SciPy 的旧版本兼容,返回值表现得像一个长度为 4 的 namedtuple,具有字段 slopeinterceptlow_slopehigh_slope,因此可以继续写:

slope, intercept, low_slope, high_slope = theilslopes(y, x) 

参考文献

[1] (1,2,3)

P.K. Sen, “基于 Kendall's tau 的回归系数估计”, J. Am. Stat. Assoc., Vol. 63, pp. 1379-1389, 1968.

[2]

H. Theil, “一种秩不变的线性和多项式回归分析方法 I, II 和 III”, Nederl. Akad. Wetensch., Proc. 53:, pp. 386-392, pp. 521-525, pp. 1397-1412, 1950.

[3]

W.L. Conover, “实用非参数统计”, 第 2 版, John Wiley and Sons, 纽约, pp. 493.

[4]

zh.wikipedia.org/wiki/Theil%E2%80%93Sen%E5%9B%9E%E5%BD%92

示例

>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt 
>>> x = np.linspace(-5, 5, num=150)
>>> y = x + np.random.normal(size=x.size)
>>> y[11:15] += 10  # add outliers
>>> y[-5:] -= 7 

计算斜率、截距和 90%置信区间。为了比较,还使用 linregress 计算最小二乘拟合:

>>> res = stats.theilslopes(y, x, 0.90, method='separate')
>>> lsq_res = stats.linregress(x, y) 

绘制结果。Theil-Sen 回归线显示为红色,虚线红线表示斜率的置信区间(请注意,虚线红线不是回归的置信区间,因为截距的置信区间未包括在内)。绿色线显示最小二乘拟合以便比较。

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(x, y, 'b.')
>>> ax.plot(x, res[1] + res[0] * x, 'r-')
>>> ax.plot(x, res[1] + res[2] * x, 'r--')
>>> ax.plot(x, res[1] + res[3] * x, 'r--')
>>> ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'g-')
>>> plt.show() 

../../_images/scipy-stats-theilslopes-1.png

scipy.stats.page_trend_test

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.page_trend_test.html#scipy.stats.page_trend_test

scipy.stats.page_trend_test(data, ranked=False, predicted_ranks=None, method='auto')

执行 Page 的检验,衡量处理之间观察结果的趋势。

Page’s Test(也称为 Page 的(L)检验)在以下情况下很有用:

  • 至少有(n \geq 3)个处理,

  • (m \geq 2)个受试者观察每种处理,并且

  • 假设观察结果具有特定的顺序。

具体来说,该检验考虑的是零假设,即

[m_1 = m_2 = m_3 \cdots = m_n,]

其中(m_j)是在处理(j)下观察量的平均值,对立假设是

[m_1 \leq m_2 \leq m_3 \leq \cdots \leq m_n,]

其中至少有一处不等式是严格的。

正如[4]所指出的,Page 的(L)检验在趋势差异的替代假设下比 Friedman 检验具有更强的统计功效,因为 Friedman 的检验只考虑观察值的平均值差异而不考虑它们的顺序。而 Spearman 的(\rho)则考虑两个变量(例如燕子的飞行速度与它所携带的椰子的重量)的排名观察之间的相关性,Page 的(L)则关注观察(例如燕子的飞行速度)在几种不同处理(例如携带不同重量的五个椰子)中的趋势,即使在多个受试者(例如一个欧洲燕子和一个非洲燕子)重复观察的情况下也是如此。

参数:

data类似数组

一个(m \times n)数组;第(i)行第(j)列的元素是与受试者(i)和处理(j)对应的观察结果。默认情况下,假设列按预测均值递增的顺序排列。

ranked布尔值,可选

默认情况下,数据被假定为观察值而不是排名;将使用scipy.stats.rankdata沿axis=1对其进行排名。如果数据以排名形式提供,请传递参数True

predicted_ranks类似数组,可选

列均值的预测排名。如果未指定,默认假设列按预测均值递增的顺序排列,因此默认的predicted_ranks是([1, 2, \dots, n-1, n])。

method,可选

选择用于计算p-值的方法。以下选项可用。

  • ‘auto’:在合理时间内选择‘exact’和‘asymptotic’之间以获得合理精度的结果(默认)

  • ‘asymptotic’:将标准化的检验统计量与正态分布进行比较

  • ‘exact’:通过比较所有可能的排名排列(在零假设下,每个排列等可能)来计算精确的p-值

返回:

resPage 趋势检验结果

一个包含属性的对象:

statisticfloat

Page’s (L) 测试统计量。

pvaluefloat

相关 p-值

方法{‘渐近’, ‘精确’}

用于计算 p-值的方法

另见

rankdata, friedmanchisquare, spearmanr

注释

[1] 所述,“这里的 (n) ‘处理’ 也可以表示 (n) 个对象、事件、表演、人员或试验,按排名排序。” 同样,(m) ‘主体’ 也可以等同于能力分组、某种控制变量的分组、进行排名的评委或某种随机复制。

计算 (L) 统计量的过程,改编自 [1],如下:

  1. “预先用严谨的逻辑确定关于实验结果预测排序的适当假设。如果没有关于任何处理排序的合理依据,那么 (L) 检验不适用。”

  2. “与其他实验一样,确定在何种置信水平下你将拒绝零假设,即实验结果与单调假设不一致。”

  3. “将实验材料分类为具有 (n) 列(处理、排名对象、条件)和 (m) 行(主体、复制组、控制变量水平)的二向表。”

  4. “记录实验观察时,对每行进行排名”,例如 ranks = scipy.stats.rankdata(data, axis=1)

  5. “对每一列中的排名求和”,例如 colsums = np.sum(ranks, axis=0)

  6. “将每个列的排名总和乘以该列的预测排名”,例如 products = predicted_ranks * colsums

  7. “将所有这类乘积求和”,例如 L = products.sum()

[1] 进一步建议使用标准化统计量

[\chi_L² = \frac{\left[12L-3mn(n+1)²\right]²}{mn²(n²-1)(n+1)}]

“近似服从自由度为 1 的卡方分布。普通使用 (\chi²) 表相当于进行双侧一致性检验。如果需要进行单侧检验,几乎总是如此,则应将卡方表中的概率 减半。”

然而,这种标准化统计量不能区分观察值是与预测排名良好相关还是与预测排名反相关。因此,我们遵循 [2] 并计算标准化统计量

[\Lambda = \frac{L - E_0}{\sqrt{V_0}},]

其中 (E_0 = \frac{1}{4} mn(n+1)²) 和 (V_0 = \frac{1}{144} mn²(n+1)(n²-1)),“这在零假设下渐近地服从正态分布”。

p-值method='exact'是通过将L的观察值与所有(n!)^m可能的排名排列生成的L值进行比较而生成的。计算是使用[5]的递归方法执行的。

p-值未针对出现并列的情况进行调整。当存在并列时,报告的 'exact' p-值可能比真实的p-值稍大(即更保守)。然而,'asymptotic' p-值往往比 'exact' p-值小(即不那么保守)。

参考文献

[1] (1,2,3,4)

Ellis Batten Page,“多重处理的有序假设:线性等级的显著性检验”,美国统计协会杂志 58(301),第 216-230 页,1963 年。

[2] (1,2)

Markus Neuhauser,非参数统计检验:计算方法,CRC Press,第 150-152 页,2012 年。

[3] (1,2)

Statext LLC,“Page's L Trend Test - Easy Statistics”,Statext - 统计学习www.statext.com/practice/PageTrendTest03.php,访问于 2020 年 7 月 12 日。

[4]

“Page's Trend Test”,维基百科,WikimediaFoundation,en.wikipedia.org/wiki/Page%27s_trend_test,访问于 2020 年 7 月 12 日。

[5]

Robert E. Odeh,“两因素布局中 Page's L 统计量的精确分布”,统计学-模拟与计算,6(1),第 49-61 页,1977 年。

示例

我们使用来自[3]的例子:询问 10 名学生对三种教学方法 - 教程、讲座和研讨会 - 进行 1-5 分的评分,其中 1 分是最低分,5 分是最高分。我们已决定以 99%的置信水平拒绝零假设,支持我们的备择假设:研讨会将获得最高评分,而教程将获得最低评分。最初,数据已经列出,每行表示一个学生对三种方法的评分,顺序如下:教程、讲座、研讨会。

>>> table = [[3, 4, 3],
...          [2, 2, 4],
...          [3, 3, 5],
...          [1, 3, 2],
...          [2, 3, 2],
...          [2, 4, 5],
...          [1, 2, 4],
...          [3, 4, 4],
...          [2, 4, 5],
...          [1, 3, 4]] 

因为教程被假设为评分最低,教程排名对应的列应该排在第一位;研讨会被假设为评分最高,所以其列应该排在最后。由于这些列已按照预测均值递增的顺序排列,我们可以直接将表传递给page_trend_test

>>> from scipy.stats import page_trend_test
>>> res = page_trend_test(table)
>>> res
PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,
 method='exact') 

这个p-值表明,在零假设下,L统计量达到如此极端值的概率为 0.1819%。因为 0.1819%小于 1%,我们有证据拒绝零假设,支持我们的备择假设,在 99%的置信水平下。

(L) 统计量的值为 133.5. 为了手动验证这一点,我们对数据进行排名,使高分对应高排名,并通过平均排名解决并列情况:

>>> from scipy.stats import rankdata
>>> ranks = rankdata(table, axis=1)
>>> ranks
array([[1.5, 3\. , 1.5],
 [1.5, 1.5, 3\. ],
 [1.5, 1.5, 3\. ],
 [1\. , 3\. , 2\. ],
 [1.5, 3\. , 1.5],
 [1\. , 2\. , 3\. ],
 [1\. , 2\. , 3\. ],
 [1\. , 2.5, 2.5],
 [1\. , 2\. , 3\. ],
 [1\. , 2\. , 3\. ]]) 

我们在每列内添加排名,将总和乘以预测排名,然后求和。

>>> import numpy as np
>>> m, n = ranks.shape
>>> predicted_ranks = np.arange(1, n+1)
>>> L = (predicted_ranks * np.sum(ranks, axis=0)).sum()
>>> res.statistic == L
True 

如在 [3] 中所述,p 值的渐近近似是正态分布的生存函数,其在标准化检验统计量处的值:

>>> from scipy.stats import norm
>>> E0 = (m*n*(n+1)**2)/4
>>> V0 = (m*n**2*(n+1)*(n**2-1))/144
>>> Lambda = (L-E0)/np.sqrt(V0)
>>> p = norm.sf(Lambda)
>>> p
0.0012693433690751756 

这与上文由 page_trend_test 报告的 p 值不完全匹配。对于 (m \leq 12) 和 (n \leq 8),渐近分布并不准确,也不保守,因此 page_trend_test 根据表格的维度和 Page 原文中的建议选择了 method='exact'。若要覆盖 page_trend_test 的选择,请提供 method 参数。

>>> res = page_trend_test(table, method="asymptotic")
>>> res
PageTrendTestResult(statistic=133.5, pvalue=0.0012693433690751756,
 method='asymptotic') 

如果数据已经排名,我们可以传入 ranks 而不是 table 来节省计算时间。

>>> res = page_trend_test(ranks,             # ranks of data
...                       ranked=True,       # data is already ranked
...                       )
>>> res
PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,
 method='exact') 

假设原始数据的制表顺序与预测均值的顺序不同,比如讲座、研讨会、教程。

>>> table = np.asarray(table)[:, [1, 2, 0]] 

由于该表格的排列与假定的顺序不一致,我们可以重新排列表格或提供 predicted_ranks。请记住,预计讲座将排在中间位置,研讨会最高,教程最低,我们传递:

>>> res = page_trend_test(table,             # data as originally tabulated
...                       predicted_ranks=[2, 3, 1],  # our predicted order
...                       )
>>> res
PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,
 method='exact') 

scipy.stats.multiscale_graphcorr

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.multiscale_graphcorr.html#scipy.stats.multiscale_graphcorr

scipy.stats.multiscale_graphcorr(x, y, compute_distance=<function _euclidean_dist>, reps=1000, workers=1, is_twosamp=False, random_state=None)

计算多尺度图相关(MGC)检验统计量。

具体而言,对于每个点,MGC 找到一个属性的k个最近邻(例如云密度),和另一个属性的l个最近邻(例如草湿度)[1]。这对(k, l)被称为“尺度”。然而,事先不知道哪些尺度会最具信息性。因此,MGC 计算所有距离对,然后有效地计算所有尺度的距离相关性。局部相关性显示哪些尺度相对于关系是最具信息性的。因此,成功发现和解释不同数据模态之间关系的关键是自适应确定哪些尺度最具信息性,以及最具信息性尺度的几何含义。这不仅提供了是否模态相关的估计,还揭示了如何进行该决定的见解。在高维数据中尤为重要,因为简单的可视化无法揭示关系给肉眼。特别是,这一实现的表征已经从[2]中得出,并在内部进行了基准测试。

参数:

x, y ndarray

如果xy的形状为(n, p)(n, q),其中n是样本数,pq是维度数,则将运行 MGC 独立性检验。另外,如果xy的形状为(n, n),并且它们是距离或相似性矩阵,则compute_distance必须发送到None。如果xy的形状为(n, p)(m, p),则将运行不配对双样本 MGC 检验。

compute_distance可调用对象,可选

计算每个数据矩阵中样本之间的距离或相似性的函数。如果xy已经是距离矩阵,则设置为None。默认使用欧氏距离度量。如果调用自定义函数,请先创建距离矩阵或创建形如compute_distance(x)的函数,其中x是计算成对距离的数据矩阵。

reps整数,可选

使用排列测试估计零假设时的复制次数。默认为1000

workers整数或类似映射的可调用对象,可选

如果 workers 是一个整数,那么将人群细分为 workers 部分,并并行评估(使用 multiprocessing.Pool <multiprocessing>)。提供 -1 来使用所有可用于进程的核心。或者提供一个类似映射的可调用对象,例如 multiprocessing.Pool.map 用于并行评估 p 值。此评估作为 workers(func, iterable) 进行。要求 func 可以被 pickle。默认为 1

is_twosampbool, optional

如果 True,将运行双样本检验。如果 xy 的形状为 (n, p)(m, p),则此选项将被覆盖并设置为 True。如果 xy 都具有形状 (n, p),并且希望运行双样本检验,则设置为 True。默认为 False。请注意,如果输入为距离矩阵,则不会运行此操作。

random_state{None, int, numpy.random.Generator,

numpy.random.RandomState}, optional

如果 seed 为 None(或 np.random),则使用 numpy.random.RandomState 单例。如果 seed 是一个整数,则使用一个新的 RandomState 实例,其种子为 seed。如果 seed 已经是 GeneratorRandomState 实例,则使用该实例。

返回:

resMGCResult

包含属性的对象:

statisticfloat

样本 MGC 测试统计量位于 [-1, 1]

pvaluefloat

通过置换获得的 p 值。

mgc_dictdict

包含额外有用结果:

  • mgc_mapndarray

  • 关系的潜在几何的二维表示。

  • opt_scale(int, int)

  • 估计的最优尺度为 (x, y) 对。

  • null_distlist

  • 来自置换矩阵的空分布。

另请参见

pearsonr

Pearson 相关系数和用于测试非相关性的 p 值。

kendalltau

计算 Kendall's tau。

spearmanr

计算 Spearman 秩相关系数。

注释

MGC 过程及其在神经科学数据上的应用的描述可在 [1] 中找到。它通过以下步骤执行:

  1. 计算并修改为零均值列的两个距离矩阵 (D^X) 和 (D^Y)。这导致两个 (n \times n) 距离矩阵 (A) 和 (B)(中心化和无偏修改) [3]

  2. 对于所有的值 (k) 和 (l),从 (1, ..., n),

    • 对于每个属性,计算 (k) 近邻图和 (l) 近邻图。这里,(G_k (i, j)) 表示 (A) 的第 (i) 行的 (k) 个最小值,(H_l (i, j)) 表示 (B) 的第 (i) 行的 (l) 个最小值

    • 让 (\circ) 表示逐元素矩阵乘积,然后使用以下统计量对局部相关性进行求和和归一化:

[c^{kl} = \frac{\sum_{ij} A G_k B H_l} {\sqrt{\sum_{ij} A² G_k \times \sum_{ij} B² H_l}}]

  1. MGC 测试统计量是 ({ c^{kl} }) 的平滑最优局部相关性。将平滑操作表示为 (R(\cdot))(本质上将所有孤立的大相关性设置为 0,将连接的大相关性保持不变),见[3]。MGC 是,

[MGC_n (x, y) = \max_{(k, l)} R \left(c^{kl} \left( x_n, y_n \right) \right)]

由于归一化,测试统计量返回一个值在 ((-1, 1)) 之间。

返回的 p 值是使用置换检验计算的。这个过程首先通过随机置换 (y) 来估计零分布,然后计算在零分布下观察到的测试统计量至少与观察到的测试统计量一样极端的概率。

MGC 需要至少 5 个样本才能获得可靠的结果。它还可以处理高维数据集。此外,通过操纵输入数据矩阵,双样本检验问题可以简化为独立性检验问题[4]。给定大小为 (p \times n) 和 (p \times m) 的样本数据 (U) 和 (V),可以如下创建数据矩阵 (X) 和 (Y):

[X = [U | V] \in \mathcal{R}^{p \times (n + m)} Y = [0_{1 \times n} | 1_{1 \times m}] \in \mathcal{R}^{(n + m)}]

然后,MGC 统计量可以像平常一样计算。这种方法可以扩展到类似的测试,比如距离相关性[4]

1.4.0 版本中的新功能。

参考文献

[1] (1,2)

Vogelstein, J. T., Bridgeford, E. W., Wang, Q., Priebe, C. E., Maggioni, M., & Shen, C. (2019). 发现和解读不同数据模态之间的关系。《ELife》。

[2]

Panda, S., Palaniappan, S., Xiong, J., Swaminathan, A., Ramachandran, S., Bridgeford, E. W., … Vogelstein, J. T. (2019). mgcpy:一个全面的高维独立性检验 Python 包。arXiv:1907.02088

[3] (1,2)

Shen, C., Priebe, C.E., & Vogelstein, J. T. (2019). 从距离相关性到多尺度图相关性。《美国统计协会杂志》。

[4] (1,2)

Shen, C. & Vogelstein, J. T. (2018). 距离和核方法在假设检验中的精确等价性。arXiv:1806.05514

示例

>>> import numpy as np
>>> from scipy.stats import multiscale_graphcorr
>>> x = np.arange(100)
>>> y = x
>>> res = multiscale_graphcorr(x, y)
>>> res.statistic, res.pvalue
(1.0, 0.001) 

要运行一个不配对的双样本检验,

>>> x = np.arange(100)
>>> y = np.arange(79)
>>> res = multiscale_graphcorr(x, y)
>>> res.statistic, res.pvalue  
(0.033258146255703246, 0.023) 

或者,如果输入的形状相同,

>>> x = np.arange(100)
>>> y = x
>>> res = multiscale_graphcorr(x, y, is_twosamp=True)
>>> res.statistic, res.pvalue  
(-0.008021809890200488, 1.0) 

scipy.stats.chi2_contingency

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.chi2_contingency.html#scipy.stats.chi2_contingency

scipy.stats.chi2_contingency(observed, correction=True, lambda_=None)

列联表中变量独立性的卡方检验。

此函数计算列联表中观察频率独立性的卡方统计量和 p 值 [1] observed 的假设检验。基于独立性假设下的边际和计算期望频率;参见 scipy.stats.contingency.expected_freq。自由度的数量使用 numpy 函数和属性表达:

dof = observed.size - sum(observed.shape) + observed.ndim - 1 

参数:

observed:array_like

列联表。表中包含每个类别中的观察频率(即出现次数)。在二维情况下,该表通常描述为“R x C 表”。

correction:bool,可选。

如果为 True,并且自由度为 1,则对连续性应用 Yates 校正。校正的效果是将每个观察值调整 0.5 向相应的期望值。

lambda_:float 或 str,可选。

默认情况下,此测试中计算的统计量是皮尔逊卡方统计量 [2]lambda_ 允许使用来自 Cressie-Read 功率差异族的统计量 [3]。有关详细信息,请参阅 scipy.stats.power_divergence

返回:

res:Chi2ContingencyResult

一个包含属性的对象:

统计量:float

测试统计量。

p 值:float

测试的 p 值。

dof:int

自由度。

expected_freq:与observed具有相同的形状。

基于表的边际和的期望频率。

另请参见

scipy.stats.contingency.expected_freq

scipy.stats.fisher_exact

scipy.stats.chisquare

scipy.stats.power_divergence

scipy.stats.barnard_exact

scipy.stats.boschloo_exact

注意事项

这种计算的有效性的常被引用的一个准则是,只有在每个单元格中的观察频率和期望频率至少为 5 时,才应使用该测试。

这是对人口不同类别独立性的检验。当观察维度为二或更多时,该检验才有意义。将该检验应用于一维表将导致期望等于观察且卡方统计量等于 0。

由于计算中存在缺失值,此函数不处理掩码数组。

scipy.stats.chisquare 一样,此函数计算卡方统计量;此函数提供的便利性在于从给定的列联表中确定预期频率和自由度。如果这些已知,并且不需要 Yates 修正,可以使用 scipy.stats.chisquare。也就是说,如果调用:

res = chi2_contingency(obs, correction=False) 

则以下为真:

(res.statistic, res.pvalue) == stats.chisquare(obs.ravel(),
                                               f_exp=ex.ravel(),
                                               ddof=obs.size - 1 - dof) 

lambda_ 参数是在 scipy 的版本 0.13.0 中添加的。

参考文献

[1]

“列联表”,zh.wikipedia.org/wiki/%E5%88%97%E8%81%94%E8%A1%A8

[2]

“皮尔逊卡方检验”,zh.wikipedia.org/wiki/%E7%9A%AE%E5%B0%94%E9%80%8A%E5%8D%A1%E6%96%B9%E6%A3%80%E9%AA%8C

[3]

Cressie, N. 和 Read, T. R. C.,“Multinomial Goodness-of-Fit Tests”,J. Royal Stat. Soc. Series B,Vol. 46, No. 3(1984),pp. 440-464。

[4]

Berger, Jeffrey S. 等人。 “Aspirin for the Primary Prevention of Cardiovascular Events in Women and Men: A Sex-Specific Meta-analysis of Randomized Controlled Trials.” JAMA, 295(3):306-313, DOI:10.1001/jama.295.3.306, 2006。

例子

[4]中,研究了阿司匹林在预防女性和男性心血管事件中的应用。研究显著结论为:

…阿司匹林疗法通过降低女性缺血性中风的风险,从而减少心血管事件的复合风险 [...]

文章列出了各种心血管事件的研究。我们将重点放在女性的缺血性中风上。

下表总结了参与者连续多年定期服用阿司匹林或安慰剂的实验结果。记录了缺血性中风的案例:

 Aspirin   Control/Placebo
Ischemic stroke     176           230
No stroke         21035         21018 

有证据表明阿司匹林减少了缺血性中风的风险吗?我们首先提出一个零假设 (H_0):

阿司匹林的效果等同于安慰剂。

让我们通过卡方检验来评估这一假设的合理性。

>>> import numpy as np
>>> from scipy.stats import chi2_contingency
>>> table = np.array([[176, 230], [21035, 21018]])
>>> res = chi2_contingency(table)
>>> res.statistic
6.892569132546561
>>> res.pvalue
0.008655478161175739 

使用 5%的显著水平,我们将拒绝零假设,支持备择假设:“阿司匹林的效果不等同于安慰剂的效果”。因为scipy.stats.contingency.chi2_contingency执行的是双侧检验,备择假设并不指示效果的方向。我们可以使用stats.contingency.odds_ratio来支持结论,即阿司匹林减少缺血性中风的风险。

下面是进一步的示例,展示如何测试更大的列联表。

一个二路示例(2 x 3):

>>> obs = np.array([[10, 10, 20], [20, 20, 20]])
>>> res = chi2_contingency(obs)
>>> res.statistic
2.7777777777777777
>>> res.pvalue
0.24935220877729619
>>> res.dof
2
>>> res.expected_freq
array([[ 12.,  12.,  16.],
 [ 18.,  18.,  24.]]) 

使用对数似然比(即“G 检验”)而不是皮尔逊卡方统计量来进行测试。

>>> res = chi2_contingency(obs, lambda_="log-likelihood")
>>> res.statistic
2.7688587616781319
>>> res.pvalue
0.25046668010954165 

一个四路示例(2 x 2 x 2 x 2):

>>> obs = np.array(
...     [[[[12, 17],
...        [11, 16]],
...       [[11, 12],
...        [15, 16]]],
...      [[[23, 15],
...        [30, 22]],
...       [[14, 17],
...        [15, 16]]]])
>>> res = chi2_contingency(obs)
>>> res.statistic
8.7584514426741897
>>> res.pvalue
0.64417725029295503 

scipy.stats.fisher_exact

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.fisher_exact.html#scipy.stats.fisher_exact

scipy.stats.fisher_exact(table, alternative='two-sided')

在 2x2 列联表上执行 Fisher 精确检验。

零假设是观察到的表的边际必须等于这些总体的边际条件下,真实几率比是一的真实几率比,并且观察是从这些总体中抽取的。返回的统计量是几率比的无条件最大似然估计,p 值是在零假设下获得至少与实际观察到的表格一样极端的概率。与 Fisher 精确检验相关的统计量和双侧 p 值定义还有其他可能的选择,请参阅注释获取更多信息。

参数:

table整数数组

2x2 列联表。元素必须是非负整数。

alternative,可选

定义备择假设。以下选项可用(默认为‘two-sided’):

  • ‘two-sided’: 底层总体的几率比不是一

  • ‘less’: 底层总体的几率比一小

  • ‘greater’: 底层总体的几率比一大

详细信息请参阅注释。

返回:

resSignificanceResult

包含属性的对象:

统计量 float

这是先前的几率比,而不是后验估计。

p 值 float

在零假设下,获得至少与实际观察到的表格一样极端的概率。

另请参见

chi2_contingency

列联表中变量独立性的卡方检验。当表中的数字较大时,可以用作fisher_exact的替代方法。

contingency.odds_ratio

计算 2x2 列联表的几率比(样本或条件极大似然估计)。

barnard_exact

Barnard 精确检验,对于 2x2 列联表来说比 Fisher 精确检验更为强大的替代方法。

boschloo_exact

Boschloo 精确检验,对于 2x2 列联表来说比 Fisher 精确检验更为强大的替代方法。

注释

零假设和 p 值

零假设是观察下层群体的真实比率为一,且这些观察是从这些群体中随机抽样的条件下成立的:结果表的边际必须与观察表的边际相等。等价地,零假设是输入表来自超几何分布,其参数为 (如 hypergeom 中所用) M = a + b + c + d, n = a + bN = a + c,其中输入表为 [[a, b], [c, d]]。这个分布的支持区间为 max(0, N + n - M) <= x <= min(N, n),或者用输入表中的值来说是 min(0, a - d) <= x <= a + min(b, c)x 可以解释为一个 2x2 表的左上元素,因此分布中的表格形式为:

[  x           n - x     ]
[N - x    M - (n + N) + x] 

例如,如果:

table = [6  2]
        [1  4] 

那么支持区间为 2 <= x <= 7,并且分布中的表格为:

[2 6]   [3 5]   [4 4]   [5 3]   [6 2]  [7 1]
[5 0]   [4 1]   [3 2]   [2 3]   [1 4]  [0 5] 

每个表格的概率由超几何分布 hypergeom.pmf(x, M, n, N) 给出。例如,这些分别是(精确到三个有效数字):

x       2      3      4      5       6        7
p  0.0163  0.163  0.408  0.326  0.0816  0.00466 

可以用以下方式计算:

>>> import numpy as np
>>> from scipy.stats import hypergeom
>>> table = np.array([[6, 2], [1, 4]])
>>> M = table.sum()
>>> n = table[0].sum()
>>> N = table[:, 0].sum()
>>> start, end = hypergeom.support(M, n, N)
>>> hypergeom.pmf(np.arange(start, end+1), M, n, N)
array([0.01631702, 0.16317016, 0.40792541, 0.32634033, 0.08158508,
 0.004662  ]) 

双侧 p 值是,在零假设下,一个随机表的概率等于或小于输入表的概率。对于我们的示例,输入表的概率(其中 x = 6)为 0.0816。概率不超过这个值的 x 值为 2、6 和 7,因此双侧 p 值为 0.0163 + 0.0816 + 0.00466 ~= 0.10256

>>> from scipy.stats import fisher_exact
>>> res = fisher_exact(table, alternative='two-sided')
>>> res.pvalue
0.10256410256410257 

对于 alternative='greater',单侧 p 值是随机表具有 x >= a 的概率,例如在我们的示例中是 x >= 6,或 0.0816 + 0.00466 ~= 0.08626

>>> res = fisher_exact(table, alternative='greater')
>>> res.pvalue
0.08624708624708627 

这相当于在 x = 5 处计算分布的生存函数(从输入表中减去 x,因为我们想要在总和中包括 x = 6 的概率):

>>> hypergeom.sf(5, M, n, N)
0.08624708624708627 

对于 alternative='less',单侧 p 值是随机表具有 x <= a 的概率(例如我们的示例中 x <= 6),或 0.0163 + 0.163 + 0.408 + 0.326 + 0.0816 ~= 0.9949

>>> res = fisher_exact(table, alternative='less')
>>> res.pvalue
0.9953379953379957 

这相当于在 x = 6 处计算分布的累积分布函数:

>>> hypergeom.cdf(6, M, n, N)
0.9953379953379957 

比率

计算得到的比率与 R 函数 fisher.test 计算的值不同。此实现返回“样本”或“无条件”最大似然估计,而 R 中的 fisher.test 使用条件最大似然估计。要计算比率的条件最大似然估计,请使用 scipy.stats.contingency.odds_ratio.

参考文献

[1]

费舍尔,罗纳德·A,“实验设计:一位女士品茶的数学。” ISBN 978-0-486-41151-4, 1935.

[2]

“费舍尔精确检验”,zh.wikipedia.org/wiki/费舍尔精确检验

[3]

Emma V. Low 等人,“确定乙酰唑胺预防急性高山病的最低有效剂量:系统评价和荟萃分析”,BMJ,345,DOI:10.1136/bmj.e6779,2012 年。

示例

在 3 中,对乙酰唑胺预防急性高山病的有效剂量进行了研究。研究显著结论如下:

每日服用 250 mg、500 mg 和 750 mg 乙酰唑胺都能有效预防急性高山病。有可用证据表明,乙酰唑胺 250 mg 是这一适应症的最低有效剂量。

以下表格总结了实验结果,一些参与者每日服用 250 mg 乙酰唑胺,而其他参与者服用安慰剂。记录了急性高山病的发病情况:

 Acetazolamide   Control/Placebo
Acute mountain sickness            7           17
No                                15            5 

有证据表明乙酰唑胺 250 mg 能减少急性高山病的风险吗?我们首先制定一个零假设 (H_0):

使用乙酰唑胺治疗和使用安慰剂的急性高山病发病几率相同。

让我们用费舍尔检验评估这一假设的可信度。

>>> from scipy.stats import fisher_exact
>>> res = fisher_exact([[7, 17], [15, 5]], alternative='less')
>>> res.statistic
0.13725490196078433
>>> res.pvalue
0.0028841933752349743 

使用 5%的显著水平,我们会拒绝零假设,支持备择假设:“与安慰剂相比,使用乙酰唑胺治疗的急性高山病发病几率较低。”

注意

因为费舍尔精确检验的零分布假设是在假定行和列的总和都是固定的情况下形成的,所以在行总和不固定的实验中应用时,其结果是保守的。

在这种情况下,列的总和是固定的;每组有 22 名受试者。但是急性高山病的发病例数却不是(也不能在进行实验前被固定)。这是一个结果。

博斯洛检验不依赖于行总和固定的假设,因此在这种情况下提供了更强大的检验。

>>> from scipy.stats import boschloo_exact
>>> res = boschloo_exact([[7, 17], [15, 5]], alternative='less')
>>> res.statistic
0.0028841933752349743
>>> res.pvalue
0.0015141406667567101 

我们验证 p 值小于fisher_exact

scipy.stats.barnard_exact

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.barnard_exact.html#scipy.stats.barnard_exact

scipy.stats.barnard_exact(table, alternative='two-sided', pooled=True, n=32)

对一个 2x2 列联表执行 Barnard 精确检验。

参数:

table 整数的 array_like

一个 2x2 列联表。元素应为非负整数。

alternative {'two-sided', 'less', 'greater'},可选

定义零假设和备择假设。默认为“双侧”。请参阅下面注释部分中的解释。

pooled 布尔值,可选

是否计算具有混合方差(如学生 t 检验中)或非混合方差(如韦尔奇 t 检验中)的分数统计。默认为True

n 整数,可选

用于构建采样方法的采样点数。请注意,由于使用scipy.stats.qmc.Sobol选择样本点,此参数将自动转换为下一个更高的 2 次幂。默认值为 32。必须为正。在大多数情况下,32 个点足以达到良好的精度。更多的点会带来性能成本。

返回:

ber BarnardExactResult

一个结果对象,具有以下属性。

统计值 浮点数

与用户选择的pooled相对应的具有混合或非混合方差的 Wald 统计量。

p 值浮点数

P 值,即在假设原假设为真的情况下,获得至少与实际观察到的分布一样极端的概率。

另请参见

chi2_contingency

列联表中变量独立性的卡方检验。

fisher_exact

一个 2x2 列联表的 Fisher 精确检验。

boschloo_exact

Boschloo 的 2x2 列联表的精确检验,这是比 Fisher 精确检验更强大的替代方法。

注释

Barnard 检验是用于分析列联表的精确检验。它检验两个分类变量的关联,并且对于 2x2 列联表而言,比 Fisher 精确检验更具有力量。

让我们定义 (X_0) 为一个 2x2 矩阵,表示观察样本,其中每列存储二项实验,如下例所示。我们还定义 (p_1, p_2) 为 (x_{11}) 和 (x_{12}) 的理论二项概率。当使用 Barnard 精确检验时,我们可以断言三种不同的零假设:

  • (H_0 : p_1 \geq p_2) 对 (H_1 : p_1 < p_2),其中 alternative = “less”

  • (H_0 : p_1 \leq p_2) 对 (H_1 : p_1 > p_2),其中 alternative = “greater”

  • (H_0 : p_1 = p_2) 对 (H_1 : p_1 \neq p_2),其中 alternative = “two-sided”(默认值)

为了计算 Barnard's 精确检验,我们使用带有汇总或非汇总方差的 Wald 统计量 [3]。在默认假设下,即两个方差相等(pooled = True),统计量计算如下:

[T(X) = \frac{ \hat{p}_1 - \hat{p}_2 }{ \sqrt{ \hat{p}(1 - \hat{p}) (\frac{1}{c_1} + \frac{1}{c_2}) } }]

其中(\hat{p}_1, \hat{p}_2)和(\hat{p})分别是(p_1, p_2)和(p)的估计量,后者是联合概率,假设(p_1 = p_2)。

如果这个假设无效(pooled = False),则统计量为:

[T(X) = \frac{ \hat{p}_1 - \hat{p}_2 }{ \sqrt{ \frac{\hat{p}_1 (1 - \hat{p}_1)}{c_1} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{c_2} } }]

然后计算 p 值:

[\sum \binom{c_1}{x_{11}} \binom{c_2}{x_{12}} \pi^{x_{11} + x_{12}} (1 - \pi)^{t - x_{11} - x_{12}}]

在此和所有 2x2 列联表(X)的和上,如下:alternative* = "less"时,* T(X) \leq T(X_0) alternative = "greater"时,* T(X) \geq T(X_0) ,或者 * T(X) \geq |T(X_0)| * 当alternative* = "two-sided"。上面,(c_1, c_2)是第 1 和 2 列的和,(t)是总和(4 个样本元素的和)。

返回的 p 值是在烦扰参数(\pi)上取的最大 p 值,其中(0 \leq \pi \leq 1)。

此函数的复杂度为(O(n c_1 c_2)),其中n是样本点的数量。

参考文献

[1]

Barnard, G. A. “2x2 表的显著性检验”。 Biometrika。 34.1/2 (1947): 123-138. DOI:dpgkg3

[2] (1,2)

Mehta, Cyrus R., 和 Pralay Senchaudhuri. “比较两个二项分布的条件与非条件精确检验”。 Cytel Software Corporation 675 (2003): 1-5.

[3]

“Wald 检验”。 维基百科en.wikipedia.org/wiki/Wald_test

例子

[2]中展示了 Barnard's 检验的一个示例。

考虑疫苗有效性研究的以下示例(Chan, 1998)。在一个 30 名受试者的随机临床试验中,15 名接种了重组 DNA 流感疫苗,另外 15 名接种了安慰剂。安慰剂组中的 15 名受试者中有 12 名最终感染了流感,而对于疫苗组,只有 15 名受试者中的 7 名(47%)感染了流感。数据表现为一个 2 x 2 表格:

 Vaccine  Placebo
Yes     7        12
No      8        3 

在进行统计假设检验时,通常使用阈值概率或显著水平来决定是否拒绝零假设(H_0)。假设我们选择了常见的显著性水平 5%。

我们的备择假设是,疫苗将降低感染该病毒的概率;即,接种疫苗后感染病毒的概率(p_1)将小于未接种疫苗后感染病毒的概率(p_2)。因此,我们使用barnard_exact选项alternative="less"调用:

>>> import scipy.stats as stats
>>> res = stats.barnard_exact([[7, 12], [8, 3]], alternative="less")
>>> res.statistic
-1.894...
>>> res.pvalue
0.03407... 

在零假设下,即疫苗不会降低感染几率的情况下,获得至少与观察数据一样极端的测试结果的概率约为 3.4%。由于这个 p 值小于我们选择的显著性水平,我们有证据来拒绝 (H_0),支持备择假设。

假设我们使用了费舍尔精确检验:

>>> _, pvalue = stats.fisher_exact([[7, 12], [8, 3]], alternative="less")
>>> pvalue
0.0640... 

在相同的显著性阈值 5%下,我们无法拒绝零假设,支持备择假设。正如在[2]中所述,巴纳德检验比费舍尔精确检验更具统计功效,因为巴纳德检验不依赖于任何边际条件。费舍尔检验应仅在两组边际都固定的情况下使用。

scipy.stats.boschloo_exact

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.boschloo_exact.html#scipy.stats.boschloo_exact

scipy.stats.boschloo_exact(table, alternative='two-sided', n=32)

在二维列联表上执行 Boschloo 的精确检验。

参数:

table整数数组

一个二维列联表。元素应为非负整数。

alternative,可选

定义了零假设和备择假设。默认为 'two-sided'。请参见下面的注释部分中的解释。

nint,可选

用于构建抽样方法中使用的抽样点数量。请注意,由于使用 scipy.stats.qmc.Sobol 选择抽样点,因此此参数将自动转换为更高的 2 的幂次方。默认为 32。必须为正数。在大多数情况下,32 个点足以达到良好的精度。更多点会带来性能成本。

返回:

berBoschlooExactResult

一个结果对象,具有以下属性。

统计量:浮点数

Boschloo 检验中使用的统计量;即来自 Fisher 精确检验的 P 值。

P 值:浮点数

P 值,即在假设原假设成立的情况下,观察到至少与实际观察到的分布一样极端的分布的概率。

参见

chi2_contingency

二维列联表中变量独立性的卡方检验。

fisher_exact

Fisher 精确检验在二维列联表上的应用。

barnard_exact

Barnard 的精确检验,它是二维列联表中比 Fisher 精确检验更强大的替代方法。

注释

Boschloo 的检验是用于分析列联表的精确检验。它检验两个分类变量之间的关联,并且是二维列联表中比 Fisher 精确检验更具统一更强大的替代方法。

Boschloo 的精确检验使用 Fisher 精确检验的 P 值作为统计量,而 Boschloo 的 P 值是在零假设下观察到这种统计量的极端值的概率。

让我们定义 (X_0) 为一个表示观察样本的二维矩阵,其中每列存储二项式实验,如下例所示。让我们还定义 (p_1, p_2) 为 (x_{11}) 和 (x_{12}) 的理论二项式概率。在使用 Boschloo 精确检验时,我们可以提出三种不同的备择假设:

  • (H_0 : p_1=p_2) 对 (H_1 : p_1 < p_2),alternative = “less”

  • (H_0 : p_1=p_2) 对 (H_1 : p_1 > p_2),alternative = “greater”

  • (H_0 : p_1=p_2) 对 (H_1 : p_1 \neq p_2),alternative = “two-sided”(默认)

当空值分布不对称时,计算双边 p 值的多种约定。在这里,我们应用这样一种约定,即双边检验的 p 值是单边检验 p 值的两倍(截断为 1.0)。请注意,fisher_exact遵循不同的约定,因此对于给定的tableboschloo_exact报告的统计量可能与fisher_exact报告的 p 值不同,当alternative='two-sided'时。

新版本 1.7.0 中的新增内容。

参考文献

[1]

R.D. Boschloo,“在检验两个概率相等时提升 2 x 2 表的条件显著水平”,Statistica Neerlandica,24(1),1970 年

[2]

“Boschloo's test”,维基百科,en.wikipedia.org/wiki/Boschloo%27s_test

[3]

Lise M. Saari 等人,“员工态度和工作满意度”,人力资源管理,43(4),395-407,2004 年,DOI:10.1002/hrm.20032

示例

在下面的例子中,我们考虑了文章“员工态度和工作满意度”[3],该文章报告了对 63 名科学家和 117 名大学教授进行的调查结果。在 63 名科学家中,有 31 名表示他们对工作非常满意,而在 117 名大学教授中,有 74 名表示他们对工作非常满意。这是否是大学教授比科学家更满意他们的工作的重要证据?下表总结了上述数据:

 college professors   scientists
Very Satisfied   74                     31
Dissatisfied     43                     32 

在进行统计假设检验时,我们通常会选择一个阈值概率或显著性水平,用来决定是否拒绝零假设(H_0)。假设我们选择常见的 5%显著性水平。

我们的备择假设是大学教授对他们的工作更满意,而不是科学家。因此,我们期望(p_1)非常满意的大学教授的比例要大于(p_2),即非常满意的科学家的比例。因此,我们调用boschloo_exact并选择alternative="greater"选项:

>>> import scipy.stats as stats
>>> res = stats.boschloo_exact([[74, 31], [43, 32]], alternative="greater")
>>> res.statistic
0.0483...
>>> res.pvalue
0.0355... 

在零假设下,即科学家比大学教授在工作中更快乐,获得至少与观察数据一样极端测试结果的概率约为 3.55%。由于此 p 值小于我们选择的显著性水平,我们有证据拒绝(H_0),支持备择假设。

scipy.stats.ttest_ind_from_stats

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ttest_ind_from_stats.html#scipy.stats.ttest_ind_from_stats

scipy.stats.ttest_ind_from_stats(mean1, std1, nobs1, mean2, std2, nobs2, equal_var=True, alternative='two-sided')

两个独立样本的均值 t 检验,从描述统计学数据。

这是检验两个独立样本具有相同平均(期望)值的零假设的检验。

参数:

mean1array_like

样本 1 的均值。

std1array_like

样本 1 的修正样本标准差(即ddof=1)。

nobs1array_like

样本 1 的观察次数。

mean2array_like

样本 2 的均值。

std2array_like

样本 2 的修正样本标准差(即ddof=1)。

nobs2array_like

样本 2 的观察次数。

equal_varbool, optional

如果为 True(默认),执行假设总体方差相等的标准独立两样本检验[1]。如果为 False,执行不假设总体方差相等的 Welch's t 检验[2]

alternative, optional

定义备择假设。以下选项可用(默认为‘two-sided’):

  • ‘two-sided’: 分布的均值不相等。

  • ‘less’: 第一个分布的平均值小于第二个分布的平均值。

  • ‘greater’: 第一个分布的平均值大于第二个分布的平均值。

1.6.0 版中的新功能。

返回:

statisticfloat or array

计算得到的 t 统计量。

pvaluefloat or array

双侧 p 值。

另请参阅

scipy.stats.ttest_ind

注释

统计量计算为(mean1 - mean2)/se,其中se为标准误差。因此,当mean1大于mean2时,统计量为正;当mean1小于mean2时,统计量为负。

此方法不会检查std1std2的任何元素是否为负数。如果在调用此方法时std1std2的任何元素为负数,则此方法将返回与分别传递numpy.abs(std1)numpy.abs(std2)相同的结果;不会抛出异常或警告。

参考文献

[1]

[zh.wikipedia.org/wiki/T 檢定#獨立樣本 t 檢定](https://zh.wikipedia.org/wiki/T 檢定#獨立樣本 t 檢定)

[2]

[zh.wikipedia.org/wiki/Welch's_t 检验](https://zh.wikipedia.org/wiki/Welch's_t 检验)

示例

假设我们有两个样本的汇总数据如下(其中样本方差为修正的样本方差):

 Sample   Sample
           Size   Mean   Variance
Sample 1    13    15.0     87.5
Sample 2    11    12.0     39.0 

对这些数据应用 t 检验(假设总体方差相等):

>>> import numpy as np
>>> from scipy.stats import ttest_ind_from_stats
>>> ttest_ind_from_stats(mean1=15.0, std1=np.sqrt(87.5), nobs1=13,
...                      mean2=12.0, std2=np.sqrt(39.0), nobs2=11)
Ttest_indResult(statistic=0.9051358093310269, pvalue=0.3751996797581487) 

对比起来,这是摘要统计数据来自的数据。利用这些数据,我们可以使用scipy.stats.ttest_ind计算相同的结果:

>>> a = np.array([1, 3, 4, 6, 11, 13, 15, 19, 22, 24, 25, 26, 26])
>>> b = np.array([2, 4, 6, 9, 11, 13, 14, 15, 18, 19, 21])
>>> from scipy.stats import ttest_ind
>>> ttest_ind(a, b)
Ttest_indResult(statistic=0.905135809331027, pvalue=0.3751996797581486) 

假设我们有二进制数据,并希望应用 t 检验来比较两个独立组中 1 的比例:

 Number of    Sample     Sample
            Size    ones        Mean     Variance
Sample 1    150      30         0.2        0.161073
Sample 2    200      45         0.225      0.175251 

样本均值 (\hat{p}) 是样本中 1 的比例,而二进制观察的方差由 (\hat{p}(1-\hat{p})) 估算。

>>> ttest_ind_from_stats(mean1=0.2, std1=np.sqrt(0.161073), nobs1=150,
...                      mean2=0.225, std2=np.sqrt(0.175251), nobs2=200)
Ttest_indResult(statistic=-0.5627187905196761, pvalue=0.5739887114209541) 

对比起来,我们可以使用 0 和 1 的数组以及scipy.stat.ttest_ind计算 t 统计量和 p 值,就像上面一样。

>>> group1 = np.array([1]*30 + [0]*(150-30))
>>> group2 = np.array([1]*45 + [0]*(200-45))
>>> ttest_ind(group1, group2)
Ttest_indResult(statistic=-0.5627179589855622, pvalue=0.573989277115258) 

scipy.stats.poisson_means_test

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.poisson_means_test.html#scipy.stats.poisson_means_test

scipy.stats.poisson_means_test(k1, n1, k2, n2, *, diff=0, alternative='two-sided')

执行泊松均值检验,又称“E-测试”。

This is a test of the null hypothesis that the difference between means of two Poisson distributions is diff. The samples are provided as the number of events k1 and k2 observed within measurement intervals (e.g. of time, space, number of observations) of sizes n1 and n2.

Parameters:

k1int

Number of events observed from distribution 1.

n1: float

Size of sample from distribution 1.

k2int

Number of events observed from distribution 2.

n2float

Size of sample from distribution 2.

difffloat, default=0

The hypothesized difference in means between the distributions underlying the samples.

alternative, optional

Defines the alternative hypothesis. The following options are available (default is ‘two-sided’):

  • ‘two-sided’: the difference between distribution means is not equal to diff
  • ‘less’: the difference between distribution means is less than diff
  • ‘greater’: the difference between distribution means is greater than diff

Returns:

statisticfloat

测试统计量(见[1] 方程式 3.3)。

pvaluefloat

The probability of achieving such an extreme value of the test statistic under the null hypothesis.

Notes

Let:

[X_1 \sim \mbox{Poisson}(\mathtt{n1}\lambda_1)]

be a random variable independent of

[X_2 \sim \mbox{Poisson}(\mathtt{n2}\lambda_2)]

and let k1 and k2 be the observed values of (X_1) and (X_2), respectively. Then poisson_means_test uses the number of observed events k1 and k2 from samples of size n1 and n2, respectively, to test the null hypothesis that

[H_0: \lambda_1 - \lambda_2 = \mathtt{diff}]

A benefit of the E-test is that it has good power for small sample sizes, which can reduce sampling costs [1]. It has been evaluated and determined to be more powerful than the comparable C-test, sometimes referred to as the Poisson exact test.

References

[1] (1,2)

Krishnamoorthy, K., & Thomson, J. (2004). A more powerful test for comparing two Poisson means. Journal of Statistical Planning and Inference, 119(1), 23-35.

[2]

Przyborowski, J., & Wilenski, H. (1940). Homogeneity of results in testing samples from Poisson series: With an application to testing clover seed for dodder. Biometrika, 31(3/4), 313-323.

Examples

假设一个园艺师希望测试从种子公司购买的苜蓿种子袋中的病草(杂草)种子数量。先前已经确定苜蓿中病草种子的数量服从泊松分布。

从袋子中取出 100 克样本,并在运送给园丁之前进行分析。样本经分析后发现不含有爬根藤种子;也就是说,k1为 0。然而,园丁到货后又从袋中取出 100 克样本。这次,在样本中发现了三颗爬根藤种子;也就是说,k2为 3。园丁想知道这种差异是否显著且不是由于偶然因素引起的。零假设是两个样本之间的差异仅仅是由于偶然因素引起的,即 (\lambda_1 - \lambda_2 = \mathtt{diff}),其中 (\mathtt{diff} = 0)。备择假设是差异不是由偶然因素引起的,即 (\lambda_1 - \lambda_2 \ne 0)。园丁选择了 5%的显著水平,以拒绝零假设,支持备择假设[2]

>>> import scipy.stats as stats
>>> res = stats.poisson_means_test(0, 100, 3, 100)
>>> res.statistic, res.pvalue
(-1.7320508075688772, 0.08837900929018157) 

P 值为 0.088,表明在零假设下观察到测试统计量的值的几率接近 9%。这超过了 5%,因此园丁不拒绝零假设,因为在这个水平上不能认为差异是显著的。

scipy.stats.ttest_ind

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ttest_ind.html#scipy.stats.ttest_ind

scipy.stats.ttest_ind(a, b, axis=0, equal_var=True, nan_policy='propagate', permutations=None, random_state=None, alternative='two-sided', trim=0, *, keepdims=False)

计算两个独立样本的均值的 T 检验。

这是一个测试,用于检验两个独立样本的均值(期望值)是否相同的空假设。该测试默认假设总体具有相同的方差。

参数:

a, b:数组类型

数组必须具有相同的形状,除了与 axis 对应的维度(默认为第一维)。

axis:整数或 None,默认为 0

如果是整数,则计算输入的轴(例如行)上的统计量。输入的每个轴切片的统计量将出现在输出的相应元素中。如果为 None,则在计算统计量之前会对输入进行展平。

equal_var:布尔值,可选

如果为 True(默认),执行假设两个独立样本具有相等总体方差的标准独立 2 样本测试 [1]。如果为 False,则执行威尔奇 t 检验,该检验不假设相等的总体方差 [2]

自版本 0.11.0 新增。

nan_policy:{‘propagate’, ‘omit’, ‘raise’}

定义如何处理输入中的 NaN。

  • propagate:如果在计算统计量的轴切片(例如行)中存在 NaN,则输出的相应条目将为 NaN。

  • omit:在执行计算时将 NaN 剔除。如果在计算统计量的轴切片中剩余的数据不足,则输出的相应条目将为 NaN。

  • raise:如果存在 NaN,则会引发 ValueError

permutations:非负整数、np.inf 或 None(默认),可选

如果为 0 或 None(默认),则使用 t 分布计算 p 值。否则,permutations 是用于使用排列检验估计 p 值的随机排列次数。如果 permutations 等于或超过池化数据的不同分区数,则执行精确测试(即每个不同分区仅使用一次)。有关详细信息,请参阅注释。

自版本 1.7.0 新增。

random_state:{None, 整数, numpy.random.Generator

numpy.random.RandomState,可选

如果 seed 为 None(或 np.random),则使用 numpy.random.RandomState 单例。如果 seed 是整数,则使用新的 RandomState 实例,并使用 seed 进行种子化。如果 seed 已经是 GeneratorRandomState 实例,则使用该实例。

用于生成排列的伪随机数生成器状态(仅在 permutations 不为 None 时使用)。

1.7.0 版本中的新功能。

alternative,可选

定义了备择假设。以下选项可用(默认为‘双侧’):

  • ‘two-sided’:样本底层分布的均值不相等。

  • ‘less’:第一个样本潜在分布的均值小于第二个样本潜在分布的均值。

  • ‘greater’:第一个样本潜在分布的平均值大于第二个样本潜在分布的平均值。

1.6.0 版本中的新功能。

trimfloat,可选

如果非零,执行修剪(Yuen’s)t 检验。定义从输入样本的每端修剪的元素比例。如果为 0(默认),则不会从任何一侧修剪元素。每个尾部修剪元素的数量是修剪次数乘以元素数量的地板值。有效范围为 0,.5)。

1.7 版本中的新功能。

keepdimsbool,默认值:False

如果设置为 True,则减少的轴作为维度大小为一的结果保留在结果中。使用此选项,结果将正确广播到输入数组。

返回:

result[TtestResult

具有以下属性的对象:

statisticfloat 或 ndarray

t 统计量。

pvaluefloat 或 ndarray

与给定备择假设相关联的 p 值。

dffloat 或 ndarray

用于计算 t 统计量的自由度数。对于排列 t 检验,此值始终为 NaN。

1.11.0 版本中的新功能。

该对象还具有以下方法:

confidence_interval(confidence_level=0.95)

为给定置信水平计算两总体均值差异的置信区间。置信区间以namedtuple返回,具有lowhigh字段。执行排列 t 检验时,不计算置信区间,lowhigh字段包含 NaN。

1.11.0 版本中的新功能。

注释

假设我们观察到两个独立样本,例如花瓣长度,并且我们正在考虑这两个样本是否来自同一总体(例如同一种花或具有相似花瓣特征的两种物种)或两个不同的总体。

t 检验量化两样本算术均值之间的差异。p 值量化在假设空白假设为真的情况下观察到的或更极端值的概率,即样本来自具有相同总体均值的总体。大于所选阈值的 p 值(例如 5%或 1%)表示我们的观察不太可能是偶然发生的。因此,我们不拒绝等同总体均值的空白假设。如果 p 值小于我们的阈值,则我们有反对等同总体均值的空白假设的证据。

默认情况下,通过将观察数据的 t 统计量与理论 t 分布进行比较来确定 p 值。当 1 < permutations < binom(n, k) 时,其中

  • ka 中的观察次数,

  • nab 中的总观察数,

  • binom(n, k) 是二项式系数(nk),

数据被汇总(连接起来),随机分配到 a 组或 b 组,并计算 t 统计量。这个过程重复进行(permutation 次),生成零假设下 t 统计量的分布,将观察数据的 t 统计量与此分布进行比较,以确定 p 值。具体来说,报告的 p 值是在 [3] 第 4.4 节中定义的 “实现显著性水平”(ASL)。请注意,还有其他使用随机置换测试估计 p 值的方法;有关其他选项,请参见更一般的 permutation_test

permutations >= binom(n, k) 时,进行精确检验:数据按每种不同方式精确分组一次。

置换检验可能计算成本高,并且不一定比分析检验更准确,但它不对基础分布的形状做出强烈假设。

常见的修剪方法被称为修剪 t 检验。有时被称为尤恩 t 检验,这是 Welch t 检验的扩展,区别在于在方差计算中使用修剪平均数以及在统计量计算中使用修剪样本大小。如果基础分布呈长尾分布或受离群值污染,建议使用修剪方法 [4]

统计量计算为 (np.mean(a) - np.mean(b))/se,其中 se 是标准误差。因此,当 a 的样本均值大于 b 的样本均值时,统计量为正;当 a 的样本均值小于 b 的样本均值时,统计量为负。

从 SciPy 1.9 开始,将不建议使用的 np.matrix 输入转换为 np.ndarray 后再执行计算。在这种情况下,输出将是一个适当形状的标量或 np.ndarray,而不是 2D 的 np.matrix。类似地,虽然忽略掩码数组的掩码元素,但输出将是一个标量或 np.ndarray,而不是具有 mask=False 的掩码数组。

参考文献

[1]

en.wikipedia.org/wiki/T-test#Independent_two-sample_t-test

[2]

en.wikipedia.org/wiki/Welch%27s_t-test

[3]

  1. Efron 和 T. Hastie. Computer Age Statistical Inference. (2016).

[4]

Yuen, Karen K. “不等总体方差的两样本修剪 t。”Biometrika,vol. 61,no. 1,1974 年,pp. 165-170。JSTOR,www.jstor.org/stable/2334299。访问日期:2021 年 3 月 30 日。

[5]

Yuen, Karen K.和 W.J. Dixon. “两样本修剪 t 的近似行为和性能。”Biometrika,vol. 60,no. 2,1973 年,pp. 369-374。JSTOR,www.jstor.org/stable/2334550。访问日期:2021 年 3 月 30 日。

示例

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng() 

具有相同均值的样本的检验:

>>> rvs1 = stats.norm.rvs(loc=5, scale=10, size=500, random_state=rng)
>>> rvs2 = stats.norm.rvs(loc=5, scale=10, size=500, random_state=rng)
>>> stats.ttest_ind(rvs1, rvs2)
Ttest_indResult(statistic=-0.4390847099199348, pvalue=0.6606952038870015)
>>> stats.ttest_ind(rvs1, rvs2, equal_var=False)
Ttest_indResult(statistic=-0.4390847099199348, pvalue=0.6606952553131064) 

ttest_ind 低估了不等方差情况下的 p 值:

>>> rvs3 = stats.norm.rvs(loc=5, scale=20, size=500, random_state=rng)
>>> stats.ttest_ind(rvs1, rvs3)
Ttest_indResult(statistic=-1.6370984482905417, pvalue=0.1019251574705033)
>>> stats.ttest_ind(rvs1, rvs3, equal_var=False)
Ttest_indResult(statistic=-1.637098448290542, pvalue=0.10202110497954867) 

n1 != n2时,等方差 t 统计量不再等于不等方差 t 统计量:

>>> rvs4 = stats.norm.rvs(loc=5, scale=20, size=100, random_state=rng)
>>> stats.ttest_ind(rvs1, rvs4)
Ttest_indResult(statistic=-1.9481646859513422, pvalue=0.05186270935842703)
>>> stats.ttest_ind(rvs1, rvs4, equal_var=False)
Ttest_indResult(statistic=-1.3146566100751664, pvalue=0.1913495266513811) 

具有不同均值、方差和 n 的 t 检验:

>>> rvs5 = stats.norm.rvs(loc=8, scale=20, size=100, random_state=rng)
>>> stats.ttest_ind(rvs1, rvs5)
Ttest_indResult(statistic=-2.8415950600298774, pvalue=0.0046418707568707885)
>>> stats.ttest_ind(rvs1, rvs5, equal_var=False)
Ttest_indResult(statistic=-1.8686598649188084, pvalue=0.06434714193919686) 

在执行排列测试时,更多的排列通常会产生更准确的结果。使用np.random.Generator来确保可重复性:

>>> stats.ttest_ind(rvs1, rvs5, permutations=10000,
...                 random_state=rng)
Ttest_indResult(statistic=-2.8415950600298774, pvalue=0.0052994700529947) 

取这两个样本,其中一个有一个极端的尾部。

>>> a = (56, 128.6, 12, 123.8, 64.34, 78, 763.3)
>>> b = (1.1, 2.9, 4.2) 

使用trim关键字执行修剪(Yuen)t 检验。例如,使用 20%修剪,trim=.2,测试将从样本a的每个尾部减少一个元素(np.floor(trim*len(a)))。它对样本b没有影响,因为np.floor(trim*len(b))为 0。

>>> stats.ttest_ind(a, b, trim=.2)
Ttest_indResult(statistic=3.4463884028073513,
 pvalue=0.01369338726499547) 

scipy.stats.mannwhitneyu

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.mannwhitneyu.html#scipy.stats.mannwhitneyu

scipy.stats.mannwhitneyu(x, y, use_continuity=True, alternative='two-sided', axis=0, method='auto', *, nan_policy='propagate', keepdims=False)

对两个独立样本执行曼-惠特尼 U 秩和检验。

曼-惠特尼 U 检验是对假设的非参数检验:样本x的底层分布与样本y的底层分布相同。它通常用作分布之间位置差异的检验。

参数:

x, y数组样本

N 维样本数组。这些数组必须是可广播的,除了给定的axis维度外。

use_continuity布尔值,可选项

是否应用连续性修正(1/2)。当method为'asymptotic'时,默认为 True;否则不起作用。

alternative,可选项

定义备择假设。默认为‘two-sided’。设F(u)G(u)xy的底层分布的累积分布函数,则可用以下备择假设:

  • ‘two-sided’: 分布不相等,即对至少一个uF(u) ≠ G(u)

  • ‘less’: x的分布在随机上小于y的分布,即F(u) > G(u)对于所有u成立。

  • ‘greater’: x的分布在随机上大于y的分布,即F(u) < G(u)对于所有u成立。

注意,上述备择假设中的数学表达式描述了底层分布的累积分布函数(CDF)。不过乍一看,不等式的方向与自然语言描述似乎不一致,但实际上并非如此。例如,假设XY是随机变量,分别服从具有 CDF FG的分布。如果对所有uF(u) > G(u),则从X抽取的样本往往小于从Y抽取的样本。

在更严格的假设集下,备择假设可以根据分布的位置来表达;见[5]第 5.1 节。

axis整数或 None,默认值:0

如果是整数,则是沿着其计算统计量的输入轴。每个轴切片(例如行)的统计量将出现在输出的相应元素中。如果为None,则在计算统计量之前将展平输入。

method,可选项

选择用于计算p-值的方法。默认为‘auto’。以下选项可用。

  • 'asymptotic': 将标准化的检验统计量与正态分布进行比较,修正并列值。

  • 'exact': 通过将观察到的U统计量与零假设下U统计量的确切分布进行比较,计算确切的p-值。不对并列值进行修正。

  • 'auto': 当一个样本的大小小于或等于 8,并且没有关系时,选择 'exact';否则选择 'asymptotic'

nan_policy

定义如何处理输入的 NaN。

  • propagate: 如果轴切片(例如行)中存在 NaN,则相应的输出条目将为 NaN。

  • omit: 在执行计算时将省略 NaN。如果沿着计算统计量的轴切片上剩余的数据不足,则相应的输出条目将为 NaN。

  • raise: 如果存在 NaN,则会引发 ValueError

keepdimsbool,默认为 False

如果设置为 True,则减少的轴将作为大小为一的维度保留在结果中。使用此选项,结果将正确地与输入数组进行广播。

返回:

resMannwhitneyuResult

包含属性的对象:

statisticfloat

与样本 x 相对应的曼-惠特尼 U 统计量。有关与样本 y 相对应的检验统计量,请参见注释。

pvaluefloat

选择的 alternative 的相关 p-value。

另请参阅

scipy.stats.wilcoxonscipy.stats.ranksumsscipy.stats.ttest_ind

注释

如果 U1 是与样本 x 对应的统计量,则与样本 y 对应的统计量为 U2 = x.shape[axis] * y.shape[axis] - U1

mannwhitneyu 用于独立样本。对于相关/配对样本,请考虑 scipy.stats.wilcoxon

method 'exact' 在没有关系且任一样本大小小于 8 时建议使用[1]。实现遵循原始提议的递归关系[1],如[3]中描述的那样。请注意,确切方法校正关系,但 mannwhitneyu 不会在数据中存在关系时引发错误或警告。

曼-惠特尼 U 检验是独立样本 t 检验的非参数版本。当来自群体的样本的均值正态分布时,请考虑 scipy.stats.ttest_ind

从 SciPy 1.9 开始,np.matrix 输入(不建议新代码使用)在执行计算之前会转换为 np.ndarray。在这种情况下,输出将是适当形状的标量或 np.ndarray,而不是 2D np.matrix。类似地,虽然忽略掩码数组的掩码元素,但输出将是标量或具有 mask=Falsenp.ndarray,而不是掩码数组。

参考文献

[1] (1,2)

H.B. Mann 和 D.R. Whitney, “On a test of whether one of two random variables is stochastically larger than the other”, 数理统计学年报, Vol. 18, pp. 50-60, 1947.

[2]

曼-惠特尼 U 检验, 维基百科, en.wikipedia.org/wiki/Mann-Whitney_U_test

[3]

A. Di Bucchianico, “Combinatorics, computer algebra, and the Wilcoxon-Mann-Whitney test”, 统计规划与推断杂志, Vol. 79, pp. 349-364, 1999.

[4] (1,2,3,4,5,6,7)

Rosie Shier, “Statistics: 2.3 The Mann-Whitney U Test”, 数学学习支持中心, 2004.

[5]

Michael P. Fay 和 Michael A. Proschan。“Wilcoxon-Mann-Whitney or t-test? On assumptions for hypothesis tests and multiple interpretations of decision rules.” 统计调查, Vol. 4, pp. 1-39, 2010. www.ncbi.nlm.nih.gov/pmc/articles/PMC2857732/

例子

我们遵循来自[4]的示例:九名随机抽样的年轻成年人被诊断为二型糖尿病的年龄如下。

>>> males = [19, 22, 16, 29, 24]
>>> females = [20, 11, 17, 12] 

我们使用曼-惠特尼 U 检验来评估男性和女性诊断年龄之间是否存在统计学显著差异。零假设是男性诊断年龄的分布与女性诊断年龄的分布相同。我们决定需要在 95%的置信水平下拒绝零假设,支持分布不同的备择假设。由于样本数非常小且数据中没有并列项,我们可以将观察到的测试统计量与零假设下测试统计量的精确分布进行比较。

>>> from scipy.stats import mannwhitneyu
>>> U1, p = mannwhitneyu(males, females, method="exact")
>>> print(U1)
17.0 

mannwhitneyu 总是报告与第一样本相关的统计量,这在本例中是男性。这与在[4]报告的 (U_M = 17) 一致。第二统计量相关的统计量可以计算:

>>> nx, ny = len(males), len(females)
>>> U2 = nx*ny - U1
>>> print(U2)
3.0 

这与在[4]报告的 (U_F = 3) 一致。双侧p-值可以从任一统计量计算,而由mannwhitneyu产生的值与在[4]报告的 (p = 0.11) 一致。

>>> print(p)
0.1111111111111111 

测试统计量的确切分布渐近正态,因此示例继续通过比较精确的p-值与使用正态近似产生的p-值。

>>> _, pnorm = mannwhitneyu(males, females, method="asymptotic")
>>> print(pnorm)
0.11134688653314041 

在这里,mannwhitneyu 报告的 p-值似乎与 [4] 给出的值 (p = 0.09) 存在冲突。原因是 [4] 没有应用由 mannwhitneyu 执行的连续性校正;mannwhitneyu 减少了测试统计量与均值 (\mu = n_x n_y / 2) 之间的距离 0.5,以校正离散统计量与连续分布的比较。在这里,使用的 (U) 统计量小于均值,因此我们在分子中加入了 0.5 以减少距离。

>>> import numpy as np
>>> from scipy.stats import norm
>>> U = min(U1, U2)
>>> N = nx + ny
>>> z = (U - nx*ny/2 + 0.5) / np.sqrt(nx*ny * (N + 1)/ 12)
>>> p = 2 * norm.cdf(z)  # use CDF to get p-value from smaller statistic
>>> print(p)
0.11134688653314041 

如果需要,我们可以禁用连续性校正,以获得与 [4] 报告的结果一致的结果。

>>> _, pnorm = mannwhitneyu(males, females, use_continuity=False,
...                         method="asymptotic")
>>> print(pnorm)
0.0864107329737 

无论我们执行精确还是渐近检验,测试统计量出现这样或更极端的概率超过 5%,因此我们不认为结果具有统计学意义。

假设在查看数据之前,我们假设女性被诊断的年龄比男性更年轻。在这种情况下,将女性年龄作为第一个输入是很自然的选择,我们将使用alternative = 'less'执行单侧检验:女性被诊断的年龄随机地小于男性。

>>> res = mannwhitneyu(females, males, alternative="less", method="exact")
>>> print(res)
MannwhitneyuResult(statistic=3.0, pvalue=0.05555555555555555) 

再次强调,在零假设下,测试统计量得到足够低的概率大于 5%,因此我们不拒绝零假设,支持我们的备择假设。

如果合理假设来自两个总体样本的均值是正态分布的,我们可以使用 t 检验进行分析。

>>> from scipy.stats import ttest_ind
>>> res = ttest_ind(females, males, alternative="less")
>>> print(res)
Ttest_indResult(statistic=-2.239334696520584, pvalue=0.030068441095757924) 

在这种假设下,p-值足够低,可以拒绝零假设,支持备择假设。

scipy.stats.bws_test

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.bws_test.html#scipy.stats.bws_test

scipy.stats.bws_test(x, y, *, alternative='two-sided', method=None)

对两个独立样本执行 Baumgartner-Weiss-Schindler 测试。

Baumgartner-Weiss-Schindler(BWS)测试是对零假设进行非参数检验,即样本x背后的分布与样本y背后的分布相同。与 Kolmogorov-Smirnov、Wilcoxon 和 Cramer-Von Mises 测试不同,BWS 测试通过差异累积分布函数(CDFs)的方差加权积分,强调分布的尾部,从而提高了许多应用中的检验能力。

参数:

x, yarray-like

1-d arrays of samples.

alternative, optional

定义备择假设。默认为‘two-sided’。设F(u)G(u)xy背后的分布的累积分布函数,则以下备择假设可用:

  • ‘two-sided’:分布不相等,即至少存在一个u使得F(u) ≠ G(u)

  • ‘less’:x背后的分布小于y背后的分布,即F(u) >= G(u)对于所有u

  • ‘greater’:x背后的分布大于y背后的分布,即F(u) <= G(u)对于所有u

在更严格的假设集下,备择假设可以用分布的位置表示;参见[2]第 5.1 节。

methodPermutationMethod, optional

配置用于计算 p 值的方法。默认是默认的PermutationMethod对象。

返回:

resPermutationTestResult

具有以下属性的对象:

statisticfloat

数据的观察检验统计量。

pvaluefloat

给定备择假设的 p 值。

null_distributionndarray

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

另见

scipy.stats.wilcoxon, scipy.stats.mannwhitneyu, scipy.stats.ttest_ind

注:

alternative=='two-sided'时,统计量由[1]第二部分中给出的方程定义。该统计量不适用于单侧备择假设;在这种情况下,统计量为[1]第二部分中给出的方程的负值。因此,当第一个样本的分布大于第二个样本的分布时,统计量趋向于为正。

参考文献

[1] (1,2,3,4,5)

Neuhäuser, M. (2005). Baumgartner-Weiss-Schindler 统计量的精确检验:一项调查。Statistical Papers, 46(1), 1-29。

[2]

Fay, M. P., & Proschan, M. A. (2010). Wilcoxon-Mann-Whitney 还是 t 检验?关于假设检验的假设和决策规则的多重解释。Statistics Surveys, 4, 1。

例如

我们遵循[1]中表 3 的示例:十四名儿童随机分为两组。他们在进行特定测试时的排名如下。

>>> import numpy as np
>>> x = [1, 2, 3, 4, 6, 7, 8]
>>> y = [5, 9, 10, 11, 12, 13, 14] 

我们使用 BWS 测试来评估两组之间是否存在统计显著差异。零假设是两组表现分布没有差异。我们决定以 1%的显著水平拒绝零假设,支持备择假设,即两组表现分布不同。由于样本量非常小,我们可以将观察到的检验统计量与在零假设下检验统计量的精确分布进行比较。

>>> from scipy.stats import bws_test
>>> res = bws_test(x, y)
>>> print(res.statistic)
5.132167152575315 

这与在[1]中报告的( B = 5.132 )一致。由bws_test产生的p-值也与在[1]中报告的( p = 0.0029 )一致。

>>> print(res.pvalue)
0.002913752913752914 

因为 p 值低于我们的 1%阈值,我们将其视为反对零假设的证据,支持备择假设,即两组表现存在差异。

scipy.stats.ranksums

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ranksums.html#scipy.stats.ranksums

scipy.stats.ranksums(x, y, alternative='two-sided', *, axis=0, nan_policy='propagate', keepdims=False)

计算两个样本的 Wilcoxon 秩和统计量。

Wilcoxon 秩和检验用于检验两组测量是否来自同一分布的零假设。备择假设是一个样本中的值更可能比另一个样本中的值大。

该检验用于比较来自连续分布的两个样本。不处理 xy 中测量之间的并列。有关并列处理和可选连续性修正,请参阅 scipy.stats.mannwhitneyu

参数:

x,yarray_like

来自两个样本的数据。

alternative, 可选

定义备择假设。默认为 ‘two-sided’。可用选项如下:

  • ‘two-sided’: xy 中的一个分布大于另一个分布。

  • ‘less’: 代表 x 分布的概率小于 y 分布的概率。

  • ‘greater’: 代表 x 分布的概率大于 y 分布的概率。

自版本 1.7.0 开始新增。

axisint 或 None,默认值为 0

如果是整数,则是输入的轴,沿着该轴计算统计量。输入的每个轴切片(例如行)的统计量将显示在输出的相应元素中。如果为 None,则在计算统计量之前将展平输入。

nan_policy

定义如何处理输入 NaN。

  • propagate: 如果在计算统计量的轴切片(例如行)中存在 NaN,则输出的相应条目将为 NaN。

  • omit: 在执行计算时将省略 NaN。如果计算统计量的轴切片中剩余的数据不足,则输出的相应条目将为 NaN。

  • raise: 如果存在 NaN,则会引发 ValueError

keepdimsbool,默认值为 False

如果设置为 True,则会将减少的轴保留在结果中,作为大小为一的维度。使用此选项,结果将正确广播到输入数组。

返回:

statisticfloat

在大样本近似下,秩和统计量被正态分布所取代的检验统计量。

pvaluefloat

检验的 p 值。

注意事项

自 SciPy 1.9 开始,np.matrix 输入(不建议新代码使用)在执行计算之前将转换为 np.ndarray。在这种情况下,输出将是适当形状的标量或 np.ndarray,而不是 2D 的 np.matrix。类似地,虽然忽略掩码数组的屏蔽元素,但输出将是标量或 np.ndarray,而不是具有 mask=False 的掩码数组。

参考文献

[1]

威尔科克森秩和检验

示例

我们可以通过计算威尔科克森秩和统计量来检验两个独立的不等大小样本是否来自同一分布。

>>> import numpy as np
>>> from scipy.stats import ranksums
>>> rng = np.random.default_rng()
>>> sample1 = rng.uniform(-1, 1, 200)
>>> sample2 = rng.uniform(-0.5, 1.5, 300) # a shifted distribution
>>> ranksums(sample1, sample2)
RanksumsResult(statistic=-7.887059,
 pvalue=3.09390448e-15) # may vary
>>> ranksums(sample1, sample2, alternative='less')
RanksumsResult(statistic=-7.750585297581713,
 pvalue=4.573497606342543e-15) # may vary
>>> ranksums(sample1, sample2, alternative='greater')
RanksumsResult(statistic=-7.750585297581713,
 pvalue=0.9999999999999954) # may vary 

p 值小于0.05表明在 5%的显著性水平下,该检验拒绝了假设。

scipy.stats.brunnermunzel

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.brunnermunzel.html#scipy.stats.brunnermunzel

scipy.stats.brunnermunzel(x, y, alternative='two-sided', distribution='t', nan_policy='propagate')

计算样本 x 和 y 上的 Brunner-Munzel 检验。

Brunner-Munzel 检验是一个非参数检验,用于检验以下原假设:当从每个组中逐个取值时,两组中获得大值的概率相等。与 Wilcoxon-Mann-Whitney U 检验不同,这不要求两组方差相同。注意,这并不假设分布相同。此检验适用于两个独立样本,可能大小不同。

参数:

x, yarray_like

样本数组,应为一维。

alternative,可选

定义备择假设。以下选项可用(默认为'双边'):

  • ‘two-sided’
  • ‘less’:单侧
  • ‘greater’:单侧

distribution,可选

定义如何获取 p 值。以下选项可用(默认为‘t’):

  • ‘t’:通过 t 分布获取 p 值
  • ‘normal’:通过标准正态分布获取 p 值。

nan_policy,可选

定义如何处理输入包含 NaN 时的情况。以下选项可用(默认为‘propagate’):

  • ‘propagate’:返回 NaN
  • ‘raise’:抛出错误
  • ‘omit’:在计算中忽略 NaN 值

返回:

statisticfloat

Brunner-Munzel 检验的 W 统计量。

pvaluefloat

假设 t 分布的 p 值。单侧或双侧,取决于 alternativedistribution 的选择。

另请参阅

mannwhitneyu

两个样本的 Mann-Whitney 秩和检验。

注释

当数据大小为 50 或更小时,Brunner 和 Munzel 建议使用 t 分布来估计 p 值。如果大小小于 10,则最好使用置换 Brunner Munzel 检验(参见[2])。

参考文献

[1]

Brunner, E. 和 Munzel, U. “非参数 Benhrens-Fisher 问题:渐近理论和小样本近似”。生物统计学期刊。Vol. 42(2000): 17-25。

[2]

Neubert, K. 和 Brunner, E. “非参数 Behrens-Fisher 问题的学生化置换检验”。计算统计与数据分析。Vol. 51(2007): 5192-5204。

示例

>>> from scipy import stats
>>> x1 = [1,2,1,1,1,1,1,1,1,1,2,4,1,1]
>>> x2 = [3,3,4,3,1,2,3,1,1,5,4]
>>> w, p_value = stats.brunnermunzel(x1, x2)
>>> w
3.1374674823029505
>>> p_value
0.0057862086661515377 

scipy.stats.mood

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.mood.html#scipy.stats.mood

scipy.stats.mood(x, y, axis=0, alternative='two-sided')

执行 Mood 的等标尺参数检验。

Mood 双样本标尺参数检验是针对两个样本是否来自具有相同标尺参数的相同分布的零假设的非参数检验。

参数:

x, yarray_like

样本数据的数组。

axisint,可选

在进行测试的轴线。xy 可以沿着 axis 有不同长度。如果 axis 是 None,则 xy 被展平,并且在展平的数组中进行测试。

alternative,可选

定义备择假设。默认为 'two-sided'。以下选项可用:

  • ‘two-sided’:xy 下 lying 分布的标尺不同。

  • ‘less’:x 下 lying 分布的标尺小于 y 下 lying 分布的标尺。

  • ‘greater’:x 下 lying 分布的标尺大于 y 下 lying 分布的标尺。

新版本 1.7.0 中引入。

返回:

resSignificanceResult

包含以下属性的对象:

statisticscalar 或 ndarray

假设检验的 z 分数。对于一维输入,返回标量。

pvaluescalar ndarray

假设检验的 p 值。

另请参阅

fligner

k 个方差相等的非参数检验

ansari

两个方差相等的非参数检验

bartlett

正态样本中 k 个方差相等的参数检验

levene

k 个方差相等的参数检验

注释

假设数据分别从概率分布 f(x)f(x/s) / s 中提取,其中 f 是某个概率密度函数。零假设是 s == 1

对于多维数组,如果输入的形状为 (n0, n1, n2, n3)(n0, m1, n2, n3),则如果 axis=1,则得到的 z 值和 p 值的形状将为 (n0, n2, n3)。注意 n1m1 不必相等,但其他维度必须相等。

参考文献

[1] Mielke, Paul W. “Note on Some Squared Rank Tests with Existing Ties.”

Technometrics,第 9 卷,第 2 期,1967 年,第 312-14 页。JSTOR,doi.org/10.2307/1266427。访问于 2022 年 5 月 18 日。

示例

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> x2 = rng.standard_normal((2, 45, 6, 7))
>>> x1 = rng.standard_normal((2, 30, 6, 7))
>>> res = stats.mood(x1, x2, axis=1)
>>> res.pvalue.shape
(2, 6, 7) 

查找标尺差异不显著的点数:

>>> (res.pvalue > 0.1).sum()
78 

以不同标尺执行测试:

>>> x1 = rng.standard_normal((2, 30))
>>> x2 = rng.standard_normal((2, 35)) * 10.0
>>> stats.mood(x1, x2, axis=1)
SignificanceResult(statistic=array([-5.76174136, -6.12650783]),
 pvalue=array([8.32505043e-09, 8.98287869e-10])) 

scipy.stats.ansari

原始文档:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ansari.html#scipy.stats.ansari

scipy.stats.ansari(x, y, alternative='two-sided', *, axis=0, nan_policy='propagate', keepdims=False)

执行 Ansari-Bradley 检验以确定尺度参数是否相等。

Ansari-Bradley 检验([1][2])是检验从两个样本抽取的分布的尺度参数相等性的非参数检验。原假设表明,x的分布的尺度与y的分布的尺度的比值为 1。

参数:

x, yarray_like

样本数据数组。

alternative,可选

定义备择假设。默认为‘two-sided’。可用的选项如下:

  • ‘two-sided’: 比例尺不等于 1。

  • ‘less’: 比例尺小于 1。

  • ‘greater’: 比例尺大于 1。

自 1.7.0 版本新增。

axisint 或 None,默认值:0

如果为整数,则输入的轴沿着计算统计量。输入的每个轴切片(例如行)的统计量将出现在输出的对应元素中。如果为None,则在计算统计量之前将对输入进行展平。

nan_policy

定义如何处理输入中的 NaN。

  • propagate: 如果轴切片(例如行)中存在 NaN,则计算统计量的对应输出将是 NaN。

  • omit: 在执行计算时将省略 NaN。如果轴切片中剩余的数据不足以进行统计计算,则对应的输出将是 NaN。

  • raise: 如果存在 NaN,则会引发ValueError

keepdimsbool,默认值:False

如果设置为 True,则会将被减少的轴保留在结果中作为大小为 1 的维度。使用此选项,结果将正确传播到输入数组。

返回:

statisticfloat

Ansari-Bradley 检验统计量。

pvaluefloat

假设检验的 p 值。

另请参阅

fligner

用于检验 k 个方差的非参数检验

mood

用于比较两个尺度参数的非参数检验

注意事项

当样本大小都小于 55 且没有任何平局时,给定的 p 值是精确的;否则,将使用 p 值的正态近似值。

自 SciPy 1.9 开始,np.matrix输入(不建议用于新代码)在执行计算之前会转换为np.ndarray。在这种情况下,输出将是适当形状的标量或np.ndarray,而不是 2D 的np.matrix。同样,虽然忽略掩码数组的掩码元素,但输出将是标量或np.ndarray,而不是具有mask=False的掩码数组。

参考文献

[1]

Ansari, A. R.和 Bradley, R. A.(1960)Dispersion 的秩和检验,数理统计学年鉴,31,1174-1189。

[2]

Sprent, Peter 和 N.C. Smeeton。应用非参数统计方法。第三版。Chapman and Hall/CRC。2001 年。第 5.8.2 节。

[3]

Nathaniel E. Helwig 的“非参数分散和平等性检验”在users.stat.umn.edu/~helwig/notes/npde-Notes.pdf

例子

>>> import numpy as np
>>> from scipy.stats import ansari
>>> rng = np.random.default_rng() 

对于这些示例,我们将创建三个随机数据集。前两个大小分别为 35 和 25,从均值为 0、标准差为 2 的正态分布中抽取。第三个数据集大小为 25,从标准差为 1.25 的正态分布中抽取。

>>> x1 = rng.normal(loc=0, scale=2, size=35)
>>> x2 = rng.normal(loc=0, scale=2, size=25)
>>> x3 = rng.normal(loc=0, scale=1.25, size=25) 

首先我们对x1x2应用ansari。这些样本来自同一分布,因此我们预计 Ansari-Bradley 检验不会导致我们得出分布比例不同的结论。

>>> ansari(x1, x2)
AnsariResult(statistic=541.0, pvalue=0.9762532927399098) 

由于 p 值接近 1,我们不能断定在比例上存在显著差异(符合预期)。

现在将测试应用于x1x3

>>> ansari(x1, x3)
AnsariResult(statistic=425.0, pvalue=0.0003087020407974518) 

在零假设相等的情况下观察到统计量极端值的概率仅为 0.03087%。我们将其视为支持备择假设的证据:从样本中抽取的分布的比例不相等。

我们可以使用alternative参数执行单侧检验。在上述示例中,x1的比例大于x3,因此x1x3的比例大于 1。这意味着当alternative='greater'时,p 值应接近 0,因此我们应该能够拒绝零假设:

>>> ansari(x1, x3, alternative='greater')
AnsariResult(statistic=425.0, pvalue=0.0001543510203987259) 

正如我们所见,p 值确实非常低。因此,使用alternative='less'应该产生较大的 p 值:

>>> ansari(x1, x3, alternative='less')
AnsariResult(statistic=425.0, pvalue=0.9998643258449039) 

scipy.stats.cramervonmises_2samp

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.cramervonmises_2samp.html#scipy.stats.cramervonmises_2samp

scipy.stats.cramervonmises_2samp(x, y, method='auto', *, axis=0, nan_policy='propagate', keepdims=False)

执行双样本 Cramér-von Mises 拟合优度检验。

这是 Cramér-von Mises 双样本检验的版本(1):对于两个独立样本 (X_1, ..., X_n) 和 (Y_1, ..., Y_m),原假设是这些样本来自相同(未指定的)连续分布。

参数:

xarray_like

1-D 数组,观测到的随机变量 (X_i) 的值。

yarray_like

1-D 数组,观测到的随机变量 (Y_i) 的值。

method,可选

用于计算 p 值的方法,请参见注意事项了解详情。默认为 ‘auto’。

axisint 或 None,默认值:0

如果是整数,则为输入的轴,沿其计算统计量。输入的每个轴切片(例如行)的统计量将显示在输出的相应元素中。如果为 None,则在计算统计量之前将对输入进行拉平。

nan_policy

定义如何处理输入的 NaN 值。

  • propagate:如果在计算统计量的轴切片(例如行)中存在 NaN,则输出的相应条目将为 NaN。

  • omit:在执行计算时,NaN 将被省略。如果在计算统计量的轴切片上剩余的数据不足,输出的相应条目将为 NaN。

  • raise:如果存在 NaN,则会引发 ValueError

keepdimsbool,默认值:False

如果设置为 True,则减少的轴将作为大小为一的维度保留在结果中。通过这个选项,结果将正确地与输入数组进行广播。

返回:

res具有属性的对象

statisticfloat

Cramér-von Mises 统计量。

pvaluefloat

p 值。

参见

cramervonmisesanderson_ksampepps_singleton_2sampks_2samp

注意事项

新版本 1.7.0 中引入。

根据 2 中的方程式 9 计算统计量。p 值的计算取决于关键字 method

  • asymptotic:通过使用检验统计量的极限分布来近似 p 值。

  • exact:通过枚举测试统计量的所有可能组合来计算精确的 p 值,参见 2。

如果 method='auto',则在两个样本包含等于或少于 20 个观测值时使用精确方法,否则使用渐近分布。

如果基础分布不是连续的,则 p 值可能是保守的(第 6.2 节在[3])。在计算检验统计量时,如果存在并列,则使用中位秩。

从 SciPy 1.9 开始,np.matrix 输入(不推荐新代码使用)在执行计算之前会转换为 np.ndarray。在这种情况下,输出将是一个相应形状的标量或 np.ndarray,而不是二维 np.matrix。类似地,虽然忽略了遮罩数组的遮罩元素,但输出将是一个标量或 np.ndarray,而不是具有 mask=False 的遮罩数组。

参考文献

[1]

en.wikipedia.org/wiki/Cramer-von_Mises_criterion

[2] (1,2)

Anderson, T.W. (1962). On the distribution of the two-sample Cramer-von-Mises criterion. The Annals of Mathematical Statistics, pp. 1148-1159.

[3]

Conover, W.J., Practical Nonparametric Statistics, 1971.

示例

假设我们希望测试由 scipy.stats.norm.rvs 生成的两个样本是否具有相同分布。我们选择显著性水平 alpha=0.05。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> x = stats.norm.rvs(size=100, random_state=rng)
>>> y = stats.norm.rvs(size=70, random_state=rng)
>>> res = stats.cramervonmises_2samp(x, y)
>>> res.statistic, res.pvalue
(0.29376470588235293, 0.1412873014573014) 

p 值超过了我们选择的显著性水平,因此我们不拒绝观察到的样本来自相同分布的原假设。

对于小样本量,可以计算精确的 p 值:

>>> x = stats.norm.rvs(size=7, random_state=rng)
>>> y = stats.t.rvs(df=2, size=6, random_state=rng)
>>> res = stats.cramervonmises_2samp(x, y, method='exact')
>>> res.statistic, res.pvalue
(0.197802197802198, 0.31643356643356646) 

基于渐近分布的 p 值是一个良好的近似,即使样本量很小。

>>> res = stats.cramervonmises_2samp(x, y, method='asymptotic')
>>> res.statistic, res.pvalue
(0.197802197802198, 0.2966041181527128) 

无论方法如何,在此示例中选择的显著性水平下,均无法拒绝原假设。

scipy.stats.epps_singleton_2samp

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.epps_singleton_2samp.html#scipy.stats.epps_singleton_2samp

scipy.stats.epps_singleton_2samp(x, y, t=(0.4, 0.8), *, axis=0, nan_policy='propagate', keepdims=False)

计算 Epps-Singleton(ES)测试统计量。

检验两个样本是否具有相同的概率分布的零假设。

参数:

x, y类似数组

要测试的两个观测样本。输入不能有多个维度。样本可以有不同的长度。

t类似数组,可选

要评估经验特征函数的点(t1, …, tn)。它应该是正的不同的数。默认值(0.4, 0.8)建议在[1]中。输入不能有多个维度。

axis整数或 None,默认:0

如果是整数,则是计算统计量时输入的轴。输入的每个轴切片(例如行)的统计量将出现在输出的相应元素中。如果None,则在计算统计量之前将对输入进行拉平处理。

nan_policy

定义如何处理输入的 NaN。

  • 传播:如果在计算统计量的轴切片(例如行)中存在 NaN,则输出的相应条目将是 NaN。

  • omit:在执行计算时将省略 NaN。如果沿计算统计量的轴切片中剩余的数据不足,则输出的相应条目将是 NaN。

  • raise:如果存在 NaN,则会引发ValueError异常。

keepdims布尔值,默认:False

如果设置为 True,则减少的轴将保留为大小为一的维度结果中。使用此选项,结果将正确广播到输入数组。

返回:

statistic浮点数

测试统计量。

pvalue浮点数

基于渐近 chi2 分布的相关 p 值。

另请参见

ks_2samp, anderson_ksamp

注意事项

在统计学中,测试两个样本是否由相同的基础分布生成是一个经典问题。广泛使用的测试是基于经验分布函数的 Kolmogorov-Smirnov(KS)测试。Epps 和 Singleton 引入了基于经验特征函数的测试[1]

与 KS 检验相比,ES 检验的一个优势是不假设连续分布。在 [1] 中,作者得出结论,该检验在许多示例中的功效也高于 KS 检验。他们建议不仅对离散样本使用 ES 检验,还建议对每个至少有 25 个观察值的连续样本使用,而对于连续情况下较小的样本量,则推荐使用 anderson_ksamp

p 值是从检验统计量的渐近分布计算得出的,该分布遵循一个 chi2 分布。如果 xy 的样本量都小于 25,那么将应用于检验统计量的小样本修正,该修正在 [1] 中提出。

t 的默认值是在 [1] 中通过考虑各种分布来确定的,并找到导致一般情况下检验功效高的良好值。在 [1] 中的表 III 给出了在该研究中测试的分布的最优值。在实现中,t 的值由半分位间距进行缩放,请参阅 [1]

从 SciPy 1.9 开始,np.matrix 输入(不建议新代码使用)在进行计算之前会被转换为 np.ndarray。在这种情况下,输出将是一个适当形状的标量或 np.ndarray,而不是二维 np.matrix。同样地,虽然忽略了掩码数组的掩码元素,但输出将是一个标量或 np.ndarray,而不是带有 mask=False 的掩码数组。

参考文献

[1] (1,2,3,4,5,6,7)

T. W. Epps 和 K. J. Singleton,“使用经验特征函数的两样本问题的综合检验”,Journal of Statistical Computation and Simulation 26,第 177–203 页,1986 年。

[2]

S. J. Goerg 和 J. Kaiser,“使用经验特征函数进行分布的非参数检验 - Epps-Singleton 两样本检验”,Stata Journal 9(3),第 454–465 页,2009 年。

scipy.stats.ks_2samp

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ks_2samp.html#scipy.stats.ks_2samp

scipy.stats.ks_2samp(data1, data2, alternative='two-sided', method='auto', *, axis=0, nan_policy='propagate', keepdims=False)

执行两样本 Kolmogorov-Smirnov 拟合优度检验。

此检验比较两个独立样本的底层连续分布 F(x) 和 G(x)。请参阅注释以了解可用的零假设和备择假设描述。

参数:

data1, data2array_like,1 维

假设两个样本观察结果数组来自连续分布,样本大小可能不同。

alternative,可选

定义零假设和备择假设。默认为 ‘two-sided’。请参阅下方注释中的解释。

method,可选

定义了计算 p 值所用的方法。以下选项可供选择(默认为 ‘auto’):

  • ‘auto’:对于小数组大小使用 ‘exact’,对于大数组使用 ‘asymp’
  • ‘exact’:使用检验统计量的确切分布
  • ‘asymp’:使用检验统计量的渐近分布

axisint 或 None,默认为 0

如果是整数,则是沿着其计算统计量的输入轴。输入的每个轴切片(例如行)的统计量将出现在输出的相应元素中。如果为 None,则在计算统计量之前会被展平。

nan_policy

定义如何处理输入的 NaN 值。

  • propagate:如果沿着计算统计量的轴切片(例如行)存在 NaN,则输出的相应条目将为 NaN。

  • omit:在执行计算时将省略 NaN。如果沿着计算统计量的轴切片中剩余的数据不足,则输出的相应条目将为 NaN。

  • raise:如果存在 NaN,则会引发 ValueError

keepdimsbool,默认为 False

如果设置为 True,则会将被减少的轴保留在结果中作为尺寸为一的维度。使用此选项,结果将正确传播至输入数组。

返回:

res:KstestResult

一个包含属性的对象:

statisticfloat

KS 检验统计量。

pvaluefloat

单尾或双尾 p 值。

statistic_locationfloat

来自 data1data2 与 KS 统计量对应的值;即,在此观察值处度量经验分布函数之间的距离。

statistic_signint

如果 data1 的经验分布函数在 statistic_location 处超过 data2 的经验分布函数,则为 +1,否则为 -1。

另请参见

kstest, ks_1samp, epps_singleton_2samp, anderson_ksamp

注意

可以使用alternative参数选择三种零假设及其对应的备择假设。

  • less: 零假设是对于所有的 x,F(x) >= G(x);备择假设是至少有一个 x 使得 F(x) < G(x)。统计量是样本的经验分布函数之间最小(最负)差异的大小。

  • greater: 零假设是对于所有的 x,F(x) <= G(x);备择假设是至少有一个 x 使得 F(x) > G(x)。统计量是样本的经验分布函数之间的最大(最正)差异。

  • two-sided: 零假设是两个分布是相同的,即对于所有的 x,F(x)=G(x);备择假设是它们不相同。统计量是样本的经验分布函数之间的最大绝对差异。

注意备择假设描述了基础分布的CDFs,而不是数据的观察值。例如,假设 x1 ~ F 和 x2 ~ G。如果对于所有的 x,F(x) > G(x),则 x1 中的值倾向于小于 x2 中的值。

如果 KS 统计量很大,则 p 值很小,这可能表明零假设被否定,支持备择假设。

如果method='exact'ks_2samp 试图计算一个精确的 p 值,即在零假设下获得与从数据计算出的测试统计值一样极端的概率。如果method='asymp',则使用渐近的 Kolmogorov-Smirnov 分布来计算近似 p 值。如果method='auto',则如果两个样本量均小于 10000,则尝试精确 p 值计算;否则,使用渐近方法。无论如何,如果尝试并失败了精确 p 值计算,将发出警告,并返回渐近 p 值。

‘two-sided’ ‘exact’ 计算计算补充概率,然后从 1 中减去。因此,它能返回的最小概率约为 1e-16。虽然算法本身是精确的,但对于大样本量,数值误差可能会累积。它最适用于其中一个样本量仅为几千的情况。

我们通常遵循 Hodges 对 Drion/Gnedenko/Korolyuk 的处理[1]

从 SciPy 1.9 开始,np.matrix 输入(不建议在新代码中使用)在执行计算之前会转换为 np.ndarray。在这种情况下,输出将是适当形状的标量或np.ndarray,而不是 2D 的np.matrix。同样地,虽然忽略掩码数组的掩码元素,但输出将是标量或np.ndarray,而不是具有mask=False的掩码数组。

参考文献

[1]

Hodges, J.L. Jr., “斯米尔诺夫双样本检验的显著性概率”,《数学档案》,3,No. 43(1958),469-486。

示例

假设我们希望检验两个样本是否来自同一分布的零假设。我们选择 95% 的置信水平;也就是说,如果 p 值小于 0.05,我们将拒绝零假设,支持备选假设。

如果第一个样本是从均匀分布抽取的,而第二个样本是从标准正态分布抽取的,我们预期将拒绝零假设。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> sample1 = stats.uniform.rvs(size=100, random_state=rng)
>>> sample2 = stats.norm.rvs(size=110, random_state=rng)
>>> stats.ks_2samp(sample1, sample2)
KstestResult(statistic=0.5454545454545454, pvalue=7.37417839555191e-15) 

实际上,p 值低于我们的阈值 0.05,因此我们拒绝零假设,支持“双边”替代假设:数据不是来自同一分布。

当两个样本来自相同分布时,我们期望数据大部分时间与零假设一致。

>>> sample1 = stats.norm.rvs(size=105, random_state=rng)
>>> sample2 = stats.norm.rvs(size=95, random_state=rng)
>>> stats.ks_2samp(sample1, sample2)
KstestResult(statistic=0.10927318295739348, pvalue=0.5438289009927495) 

正如预期的那样,p 值为 0.54 不低于我们的阈值 0.05,因此我们无法拒绝零假设。

然而,假设第一个样本是从向更大值偏移的正态分布中抽取的。在这种情况下,底层分布的累积密度函数(CDF)倾向于小于第二个样本的 CDF。因此,我们预计将拒绝零假设,采用alternative='less'

>>> sample1 = stats.norm.rvs(size=105, loc=0.5, random_state=rng)
>>> stats.ks_2samp(sample1, sample2, alternative='less')
KstestResult(statistic=0.4055137844611529, pvalue=3.5474563068855554e-08) 

而且,确实,p 值小于我们的阈值,我们拒绝零假设,支持备选假设。

scipy.stats.kstest

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.kstest.html#scipy.stats.kstest

scipy.stats.kstest(rvs, cdf, args=(), N=20, alternative='two-sided', method='auto', *, axis=0, nan_policy='propagate', keepdims=False)

执行(单样本或双样本)Kolmogorov-Smirnov 拟合优度检验。

单样本检验将样本的底层分布 F(x)与给定分布 G(x)进行比较。双样本检验比较两个独立样本的底层分布。这两个检验仅适用于连续分布。

参数:

rvs字符串、array_like 或可调用对象

如果是数组,则应该是随机变量观测的 1-D 数组。如果是可调用对象,则应该是生成随机变量的函数;它需要一个关键字参数size。如果是字符串,则应该是scipy.stats中分布的名称,将用于生成随机变量。

cdf字符串、array_like 或可调用对象

如果 array_like,则应该是随机变量观测的 1-D 数组,并执行双样本检验(rvs 必须是 array_like)。如果是可调用对象,则使用该可调用对象计算 cdf。如果是字符串,则应该是scipy.stats中分布的名称,将用作 cdf 函数。

args元组、序列,可选

分布参数,如果rvscdf是字符串或可调用对象。

N整数,可选

如果rvs为字符串或可调用对象,则为样本大小。默认值为 20。

alternative,可选

定义零假设和备择假设。默认为‘two-sided’。请参见下面的说明。

method,可选

定义用于计算 p 值的分布。提供以下选项(默认为‘auto’):

  • ‘auto’:选择其他选项之一。
  • ‘exact’:使用测试统计量的精确分布。
  • ‘approx’:用两倍的单侧概率近似计算双侧概率
  • ‘asymp’:使用测试统计量的渐近分布

axisint 或 None,默认为 0

如果是 int,则是沿着其计算统计量的输入轴(例如行)的轴。输入的每个轴切片(例如行)的统计量将显示在输出的相应元素中。如果为None,则在计算统计量之前将对输入进行拉平。

nan_policy

定义如何处理输入的 NaN。

  • propagate:如果沿着计算统计量的轴切片(例如行)存在 NaN,则输出的相应条目将为 NaN。

  • omit:执行计算时将忽略 NaN。如果沿着计算统计量的轴切片中剩余的数据不足,则输出的相应条目将为 NaN。

  • raise: 如果存在 NaN,则会引发ValueError

keepdims布尔值,默认为 False

如果设置为 True,则被减少的轴将保留在结果中作为大小为一的维度。使用此选项,结果将正确地广播到输入数组。

返回:

res:KstestResult

一个包含属性的对象:

统计量浮点数

KS 检验统计量,可以是 D+、D-或者两者中的最大值。

p 值浮点数

单侧或双侧 p 值。

statistic_location 浮点数

在单样本检验中,这是与 KS 统计量对应的rvs的值;即,在这个观察点上测量经验分布函数与假设的累积分布函数之间的距离。

在双样本检验中,这是与 KS 统计量对应的rvscdf的值;即,在这个观察值上测量经验分布函数之间的距离。

statistic_signint

在单样本检验中,如果 KS 统计量是经验分布函数与假设的累积分布函数之间的最大正差异(D+),则此值为+1;如果 KS 统计量是最大负差异(D-),则此值为-1。

在双样本检验中,如果rvs的经验分布函数在statistic_location处超过cdf的经验分布函数,则为+1;否则为-1。

另见

ks_1samp, ks_2samp

可以使用alternative参数选择三种零假设及相应的备择假设。

  • 双边检验:零假设是两个分布在所有点上相同,即 F(x)=G(x);备择假设是它们不相同。

  • 小于:零假设是对所有 x,F(x) >= G(x);备择假设是对至少一个 x,F(x) < G(x)。

  • 大于:零假设是对所有 x,F(x) <= G(x);备择假设是对至少一个 x,F(x) > G(x)。

注意备择假设描述的是底层分布的CDFs,而不是观察值。例如,假设 x1 服从 F,x2 服从 G。如果对所有 x,F(x) > G(x),则 x1 中的值倾向于小于 x2 中的值。

从 SciPy 1.9 开始,不推荐新代码使用np.matrix输入,在计算之前会被转换为np.ndarray。在这种情况下,输出将是一个适当形状的标量或np.ndarray,而不是 2D 的np.matrix。类似地,虽然被屏蔽的数组元素会被忽略,但输出将是一个标量或np.ndarray,而不是带有mask=False的屏蔽数组。

示例

假设我们希望检验样本是否按标准正态分布,我们选择 95%的置信水平;也就是说,如果 p 值小于 0.05,我们将拒绝零假设,支持备择假设。

在测试均匀分布数据时,我们预期将拒绝零假设。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> stats.kstest(stats.uniform.rvs(size=100, random_state=rng),
...              stats.norm.cdf)
KstestResult(statistic=0.5001899973268688, pvalue=1.1616392184763533e-23) 

的确,p 值低于我们的 0.05 阈值,因此我们拒绝零假设,支持默认的“双边”备择假设:数据按标准正态分布分布。

在测试来自标准正态分布的随机变量时,我们预期大部分时间数据与零假设一致。

>>> x = stats.norm.rvs(size=100, random_state=rng)
>>> stats.kstest(x, stats.norm.cdf)
KstestResult(statistic=0.05345882212970396, pvalue=0.9227159037744717) 

如预期,p 值为 0.92 不低于我们的 0.05 阈值,因此我们不能拒绝零假设。

然而,假设随机变量按向更大值偏移的正态分布分布。在这种情况下,基础分布的累积密度函数(CDF)倾向于小于标准正态分布的 CDF。因此,我们期望零假设在alternative='less'时被拒绝:

>>> x = stats.norm.rvs(size=100, loc=0.5, random_state=rng)
>>> stats.kstest(x, stats.norm.cdf, alternative='less')
KstestResult(statistic=0.17482387821055168, pvalue=0.001913921057766743) 

并且,由于 p 值小于我们的阈值,我们拒绝零假设,支持备择假设。

为了方便起见,可以使用分布名称作为第二个参数执行先前的测试。

>>> stats.kstest(x, "norm", alternative='less')
KstestResult(statistic=0.17482387821055168, pvalue=0.001913921057766743) 

上述示例都是与ks_1samp执行的单样本测试相同的。请注意kstest也可以执行与ks_2samp相同的双样本测试。例如,当两个样本来自相同分布时,我们预期大部分时间数据与零假设一致。

>>> sample1 = stats.laplace.rvs(size=105, random_state=rng)
>>> sample2 = stats.laplace.rvs(size=95, random_state=rng)
>>> stats.kstest(sample1, sample2)
KstestResult(statistic=0.11779448621553884, pvalue=0.4494256912629795) 

如预期,p 值为 0.45 不低于我们的 0.05 阈值,因此我们不能拒绝零假设。

scipy.stats.f_oneway

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.f_oneway.html#scipy.stats.f_oneway

scipy.stats.f_oneway(*samples, axis=0)

进行单因素方差分析。

单因素方差分析检验零假设是两个或多个组具有相同的总体均值。该检验适用于来自两个或多个组的样本,可能大小不同。

参数:

sample1, sample2, …array_like

每组的样本测量值。必须至少有两个参数。如果数组是多维的,则除了 axis 外,数组的所有维度必须相同。

axisint,可选

应用检验的输入数组的轴。默认为 0。

返回:

statisticfloat

测试的计算 F 统计量。

pvaluefloat

来自 F 分布的相关 p 值。

警告:

ConstantInputWarning

如果所有输入数组中的每个值都相同,则会引发错误。在这种情况下,F 统计量要么是无穷大,要么未定义,因此返回 np.infnp.nan

DegenerateDataWarning

如果任何输入数组的长度为 0,或者所有输入数组的长度为 1,则会引发错误。在这些情况下,返回 np.nan 的 F 统计量和 p 值。

注意事项

方差分析检验有重要的假设条件,这些条件必须满足才能使相关的 p 值有效。

  1. 样本是独立的。

  2. 每个样本都来自一个正态分布的总体。

  3. 所有组的总体标准差相等。这种特性称为等方差性。

如果对于给定的数据集这些假设不成立,可能仍然可以使用 Kruskal-Wallis H 检验 (scipy.stats.kruskal) 或 Alexander-Govern 测试 (scipy.stats.alexandergovern),尽管可能会有一些功效损失。

每组的长度必须至少为一,且至少有一组的长度大于一。如果这些条件不满足,会生成警告并返回 (np.nan, np.nan)。

如果每组中的所有值都相同,并且至少存在两组具有不同值,则函数会生成警告并返回 (np.inf, 0)。

如果所有组中的所有值都相同,函数会生成警告并返回 (np.nan, np.nan)。

算法来自 Heiman 的 [2],第 394-7 页。

参考文献

[1]

R. Lowry,《推断统计的概念与应用》,第十四章,2014 年,vassarstats.net/textbook/

[2]

G.W. Heiman,“理解心理学研究方法与统计学:整合导论”,霍顿·米夫林和公司,2001 年。

[3]

G.H. McDonald,“生物统计手册”,单因素方差分析。www.biostathandbook.com/onewayanova.html

例子

>>> import numpy as np
>>> from scipy.stats import f_oneway 

下面是关于贻贝 Mytilus trossulus 的壳测量数据[3](通过除以长度标准化的前附加肌瘢痕长度),来自五个地点的数据:俄勒冈州提拉穆克;俄勒冈州纽波特;阿拉斯加州彼得堡;俄罗斯马加丹;芬兰特瓦尔米内,这些数据来自 McDonald 等人(1991 年)使用的大数据集。

>>> tillamook = [0.0571, 0.0813, 0.0831, 0.0976, 0.0817, 0.0859, 0.0735,
...              0.0659, 0.0923, 0.0836]
>>> newport = [0.0873, 0.0662, 0.0672, 0.0819, 0.0749, 0.0649, 0.0835,
...            0.0725]
>>> petersburg = [0.0974, 0.1352, 0.0817, 0.1016, 0.0968, 0.1064, 0.105]
>>> magadan = [0.1033, 0.0915, 0.0781, 0.0685, 0.0677, 0.0697, 0.0764,
...            0.0689]
>>> tvarminne = [0.0703, 0.1026, 0.0956, 0.0973, 0.1039, 0.1045]
>>> f_oneway(tillamook, newport, petersburg, magadan, tvarminne)
F_onewayResult(statistic=7.121019471642447, pvalue=0.0002812242314534544) 

f_oneway 接受多维输入数组。当输入是多维的且未给定axis时,测试沿输入数组的第一个轴执行。对于以下数据,测试将执行三次,每次对应每列数据。

>>> a = np.array([[9.87, 9.03, 6.81],
...               [7.18, 8.35, 7.00],
...               [8.39, 7.58, 7.68],
...               [7.45, 6.33, 9.35],
...               [6.41, 7.10, 9.33],
...               [8.00, 8.24, 8.44]])
>>> b = np.array([[6.35, 7.30, 7.16],
...               [6.65, 6.68, 7.63],
...               [5.72, 7.73, 6.72],
...               [7.01, 9.19, 7.41],
...               [7.75, 7.87, 8.30],
...               [6.90, 7.97, 6.97]])
>>> c = np.array([[3.31, 8.77, 1.01],
...               [8.25, 3.24, 3.62],
...               [6.32, 8.81, 5.19],
...               [7.48, 8.83, 8.91],
...               [8.59, 6.01, 6.07],
...               [3.07, 9.72, 7.48]])
>>> F, p = f_oneway(a, b, c)
>>> F
array([1.75676344, 0.03701228, 3.76439349])
>>> p
array([0.20630784, 0.96375203, 0.04733157]) 

scipy.stats.tukey_hsd

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.tukey_hsd.html#scipy.stats.tukey_hsd

scipy.stats.tukey_hsd(*args)

对多个处理进行 Tukey's HSD 测试以比较均值的等性。

Tukey's HSD(Tukey 的显著差异)测试对一组样本执行均值的两两比较。而方差分析(如f_oneway)评估每个样本底层真实均值是否相同,Tukey's HSD 则是用于比较每个样本均值与其他每个样本均值的事后检验。

零假设是样本底层分布的均值相同。计算每个可能的样本配对的检验统计量,其实就是样本均值之差。对于每对,p 值是在零假设下(及其他假设;见注意事项)观察到统计量的极端值的概率,考虑到正在执行许多两两比较。还提供了每对均值差异的置信区间。

参数:

sample1, sample2, …array_like

每组的样本测量值。至少必须有两个参数。

返回:

resultTukeyHSDResult 实例

返回值是具有以下属性的对象:

statisticfloat ndarray

每次比较的测试统计量。索引(i, j)处的元素是组ij之间的统计量。

pvaluefloat ndarray

每次比较的测试 p 值。索引(i, j)处的元素是组ij之间的 p 值。

该对象具有以下方法:

confidence_interval(confidence_level=0.95):

计算指定置信水平的置信区间。

另见

dunnett

对比一组控制组的均值。

注意事项

该测试的使用依赖于几个假设。

  1. 观测值在组内和组间是独立的。

  2. 每组内和组间的观测值均服从正态分布。

  3. 从中抽取样本的分布具有相同的有限方差。

测试的原始制定是针对相等样本大小的 [6]。在样本大小不等的情况下,测试使用 Tukey-Kramer 方法 [4]

参考文献

[1]

NIST/SEMATECH 统计方法电子手册,“7.4.7.1. Tukey 方法。” www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm,2020 年 11 月 28 日。

[2]

Abdi, Herve & Williams, Lynne. (2021). “Tukey's Honestly Significant Difference (HSD) Test.” personal.utdallas.edu/~herve/abdi-HSD2010-pretty.pdf

[3]

“使用 SAS PROC ANOVA 和 PROC GLM 进行单因素方差分析.” SAS 教程, 2007, www.stattutorials.com/SAS/TUTORIAL-PROC-GLM.htm.

[4]

Kramer, Clyde Young. “扩展多重范围检验以处理具有不等复制次数的组均值.” 生物统计学, vol. 12, no. 3, 1956, pp. 307-310. JSTOR, www.jstor.org/stable/3001469. 访问于 2021 年 5 月 25 日.

[5]

NIST/SEMATECH 统计方法电子手册, “7.4.3.3. 方差分析表及均值假设检验” www.itl.nist.gov/div898/handbook/prc/section4/prc433.htm, 2021 年 6 月 2 日.

[6]

Tukey, John W. “Comparing Individual Means in the Analysis of Variance.” 生物统计学, vol. 5, no. 2, 1949, pp. 99-114. JSTOR, www.jstor.org/stable/3001913. 访问于 2021 年 6 月 14 日。

示例

这里是比较三种头痛药物的缓解时间的数据,单位为分钟。数据改编自 [3]

>>> import numpy as np
>>> from scipy.stats import tukey_hsd
>>> group0 = [24.5, 23.5, 26.4, 27.1, 29.9]
>>> group1 = [28.4, 34.2, 29.5, 32.2, 30.1]
>>> group2 = [26.1, 28.3, 24.3, 26.2, 27.8] 

我们希望查看各组均值是否显著不同。首先,通过箱线图进行视觉检查。

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(1, 1)
>>> ax.boxplot([group0, group1, group2])
>>> ax.set_xticklabels(["group0", "group1", "group2"]) 
>>> ax.set_ylabel("mean") 
>>> plt.show() 

../../_images/scipy-stats-tukey_hsd-1_00_00.png

从箱线图中,我们可以看到第 1 组到第 2 组和第 3 组的四分位数范围有重叠,但我们可以应用 tukey_hsd 测试以确定均值差异是否显著。我们设置显著水平为 .05 以拒绝零假设。

>>> res = tukey_hsd(group0, group1, group2)
>>> print(res)
Tukey's HSD Pairwise Group Comparisons (95.0% Confidence Interval)
Comparison  Statistic  p-value   Lower CI   Upper CI
(0 - 1)     -4.600      0.014     -8.249     -0.951
(0 - 2)     -0.260      0.980     -3.909      3.389
(1 - 0)      4.600      0.014      0.951      8.249
(1 - 2)      4.340      0.020      0.691      7.989
(2 - 0)      0.260      0.980     -3.389      3.909
(2 - 1)     -4.340      0.020     -7.989     -0.691 

零假设是每组具有相同的均值。对比 group0group1,以及 group1group2 的 p 值均不超过 .05,因此我们拒绝它们具有相同均值的零假设。对比 group0group2 的 p 值超过 .05,因此我们接受它们均值无显著差异的零假设。

我们还可以计算与我们选择的置信水平相关的置信区间。

>>> group0 = [24.5, 23.5, 26.4, 27.1, 29.9]
>>> group1 = [28.4, 34.2, 29.5, 32.2, 30.1]
>>> group2 = [26.1, 28.3, 24.3, 26.2, 27.8]
>>> result = tukey_hsd(group0, group1, group2)
>>> conf = res.confidence_interval(confidence_level=.99)
>>> for ((i, j), l) in np.ndenumerate(conf.low):
...     # filter out self comparisons
...     if i != j:
...         h = conf.high[i,j]
...         print(f"({i} - {j}) {l:>6.3f}  {h:>6.3f}")
(0 - 1) -9.480  0.280
(0 - 2) -5.140  4.620
(1 - 0) -0.280  9.480
(1 - 2) -0.540  9.220
(2 - 0) -4.620  5.140
(2 - 1) -9.220  0.540 
posted @ 2024-06-27 17:06  绝不原创的飞龙  阅读(13)  评论(0编辑  收藏  举报