SciPy 1.12 中文文档（六十二）

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

scipy.stats.anderson

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.anderson.html#scipy.stats.anderson

scipy.stats.anderson(x, dist='norm')

针对来自特定分布的数据的 Anderson-Darling 检验。

Anderson-Darling 检验测试零假设，即样本来自符合特定分布的总体。对于 Anderson-Darling 检验，关键值取决于正在测试的分布类型。此函数适用于正态、指数、逻辑、威布尔最小型或 Gumbel（极值型 I 型）分布。

参数：

xarray_like

样本数据数组。

dist{‘norm’, ‘expon’, ‘logistic’, ‘gumbel’, ‘gumbel_l’, ‘gumbel_r’, ‘extreme1’, ‘weibull_min’}, 可选

要测试的分布类型。默认为‘norm’。‘extreme1’、‘gumbel_l’ 和 ‘gumbel’ 是同一分布的同义词。

返回：

resultAndersonResult

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

statisticfloat

Anderson-Darling 检验统计量。

critical_valueslist

此分布的关键值。

significance_levellist

对应关键值的显著性水平，以百分比表示。函数返回针对不同分布的一组不同显著性水平的关键值。

fit_resultFitResult

包含拟合分布到数据结果的对象。

另请参见

kstest

检验拟合优度的 Kolmogorov-Smirnov 检验。

注释

提供的关键值适用于以下显著性水平：

正态/指数

15%，10%，5%，2.5%，1%

逻辑分布

25%，10%，5%，2.5%，1%，0.5%

gumbel_l / gumbel_r

25%，10%，5%，2.5%，1%

威布尔最小型

50%，25%，15%，10%，5%，2.5%，1%，0.5%

如果返回的统计量大于这些关键值，那么对应显著性水平，可以拒绝数据来自所选分布的零假设。返回的统计量在参考资料中称为“A2”。

对于weibull_min，最大似然估计已知是具有挑战性的。如果测试成功返回，则最大似然估计的一阶条件已经验证，并且临界值相对较好地对应于显著性水平，前提是样本足够大（>10 个观测值 [7]）。然而，对于一些数据，特别是没有左尾的数据，anderson可能会导致错误消息。在这种情况下，考虑使用scipy.stats.monte_carlo_test执行自定义拟合优度检验。

参考文献

[1]

www.itl.nist.gov/div898/handbook/prc/section2/prc213.htm

[2]

Stephens, M. A. (1974). 拟合优度的 EDF 统计量及其一些比较，美国统计协会杂志，第 69 卷，第 730-737 页。

[3]

Stephens, M. A. (1976). 未知参数拟合优度统计的渐近结果，统计学年鉴，第 4 卷，第 357-369 页。

[4]

Stephens, M. A. (1977). 极值分布的拟合优度，生物统计学，第 64 卷，第 583-588 页。

[5]

Stephens, M. A. (1977). 拟合优度及其与指数性测试的特别参考，技术报告编号 262，斯坦福大学统计系，斯坦福，加州。

[6]

Stephens, M. A. (1979). 基于经验分布函数的 Logistic 分布拟合优度检验，生物统计学，第 66 卷，第 591-595 页。

[7]

Richard A. Lockhart 和 Michael A. Stephens，“三参数 Weibull 分布的估计和拟合检验”，英国皇家统计学会期刊 B 系列（方法学），第 56 卷，第 3 期（1994 年），第 491-500 页，表 0。

例子

检验一个随机样本是否来自正态分布的零假设（具体均值和标准差未指定）。

>>> import numpy as np
>>> from scipy.stats import anderson
>>> rng = np.random.default_rng()
>>> data = rng.random(size=35)
>>> res = anderson(data)
>>> res.statistic
0.8398018749744764
>>> res.critical_values
array([0.527, 0.6  , 0.719, 0.839, 0.998])
>>> res.significance_level
array([15\. , 10\. ,  5\. ,  2.5,  1\. ]) 

统计量的值（勉强）超过了显著性水平为 2.5%的临界值，因此零假设可以在 2.5%的显著性水平下被拒绝，但不能在 1%的显著性水平下被拒绝。

scipy.stats.cramervonmises

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.cramervonmises.html#scipy.stats.cramervonmises

scipy.stats.cramervonmises(rvs, cdf, args=(), *, axis=0, nan_policy='propagate', keepdims=False)

执行单样本 Cramér-von Mises 拟合优度检验。

此操作用于测试累积分布函数 (F) 的拟合优度，与假定为独立同分布的观察随机变量 (X_1, ..., X_n) 的经验分布函数 (F_n) 相比较（[1]）。零假设是 (X_i) 具有累积分布 (F)。

参数：

rvsarray_like

一维数组，包含随机变量 (X_i) 的观测值。

cdfstr 或 可调用对象

用于测试观测值的累积分布函数 (F)。如果是字符串，应该是scipy.stats中分布的名称。如果是可调用对象，将使用该可调用对象来计算累积分布函数：cdf(x, *args) -> float。

args元组，可选

分布参数。假设这些是已知的；请参阅注释。

axis整数或 None，默认为 0

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

nan_policy{‘传播’, ‘省略’, ‘提高’}

定义如何处理输入的 NaN。

• 传播: 如果轴切片（例如行）中存在 NaN，则输出的相应条目将为 NaN。

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

• 提高: 如果存在 NaN，则会引发 ValueError。

keepdims布尔值，默认为 False

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

返回：

res具有属性的对象

统计量为 float

Cramér-von Mises 统计量。

p 值为 float

p 值。

参见

kstest, cramervonmises_2samp

注释

从版本 1.6.0 开始。

p 值依赖于方程式 1.8 中给出的近似值[2]。重要的是要记住，只有在测试简单假设时（即参考分布的参数已知）才能准确计算 p 值。如果参数是从数据中估计得出的（复合假设），则计算出的 p 值不可靠。

从 SciPy 1.9 开始，np.matrix 输入（不建议在新代码中使用）在执行计算之前将转换为 np.ndarray。在这种情况下，输出将是适当形状的标量或 np.ndarray，而不是二维的 np.matrix。同样地，虽然掩码数组的掩码元素被忽略，输出将是适当形状的标量或 np.ndarray，而不是具有 mask=False 的掩码数组。

参考文献

[1]

Cramér-von Mises 准则，维基百科，en.wikipedia.org/wiki/Cram%C3%A9r%E2%80%93von_Mises_criterion

[2]

Csörgő, S. 和 Faraway, J.（1996 年）。Cramér-von Mises 统计量的精确和渐近分布。《皇家统计学会杂志》，pp. 221-234。

示例

假设我们希望测试由 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=500, random_state=rng)
>>> res = stats.cramervonmises(x, 'norm')
>>> res.statistic, res.pvalue
(0.1072085112565724, 0.5508482238203407) 

P 值超过我们选择的显著性水平，因此我们不拒绝假设所观察的样本是从标准正态分布中抽取的。

现在假设我们希望检查将同样的样本移动 2.1 是否与从均值为 2 的正态分布中抽取一致。

>>> y = x + 2.1
>>> res = stats.cramervonmises(y, 'norm', args=(2,))
>>> res.statistic, res.pvalue
(0.8364446265294695, 0.00596286797008283) 

在这里，我们使用了 args 关键字来指定要对其进行数据测试的正态分布的均值（loc）。这相当于以下内容，其中我们创建一个均值为 2.1 的冻结正态分布，然后将其 cdf 方法作为参数传递。

>>> frozen_dist = stats.norm(loc=2)
>>> res = stats.cramervonmises(y, frozen_dist.cdf)
>>> res.statistic, res.pvalue
(0.8364446265294695, 0.00596286797008283) 

在任一情况下，如果 P 值小于我们选择的显著性水平，我们将拒绝假设所观察的样本是从均值为 2（默认方差为 1）的正态分布中抽取的。

scipy.stats.ks_1samp

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

scipy.stats.ks_1samp(x, cdf, args=(), alternative='two-sided', method='auto', *, axis=0, nan_policy='propagate', keepdims=False)

执行单样本 Kolmogorov-Smirnov 拟合优度检验。

此测试比较样本的基础分布 F(x) 与给定连续分布 G(x)。请参阅备注以获取可用的零假设和备择假设的描述。

参数：

x array_like

一维数组，表示 iid 随机变量的观察值。

cdf 可调用函数

用于计算 cdf 的可调用函数。

args 元组，序列，可选

与 cdf 一起使用的分布参数。

alternative {‘two-sided’, ‘less’, ‘greater’}，可选

定义零假设和备择假设。默认为 ‘two-sided’。请参阅下面的备注中的解释。

method {‘auto’, ‘exact’, ‘approx’, ‘asymp’}，可选

定义用于计算 p 值的分布。提供以下选项（默认为 ‘auto’）：

• ‘auto’：选择其他选项之一。
• ‘exact’：使用检验统计量的精确分布。
• ‘approx’：用两倍的单侧概率近似计算双侧概率
• ‘asymp’: 使用检验统计量的渐近分布

axis 整数或 None，默认为 0

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

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

定义如何处理输入的 NaN。

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

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

• raise：如果存在 NaN，则会引发 ValueError。

keepdims 布尔值，默认为 False

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

返回：

res：KstestResult

一个包含属性的对象：

statisticfloat

KS 检验统计量，可以是 D+、D- 或 D（两者中的最大值）

pvaluefloat

单尾或双尾 p 值。

statistic_locationfloat

值 x 对应于 KS 统计量；即，在此观察值处测量经验分布函数与假设的累积分布函数之间的距离。

statistic_signint

如果 KS 统计量是经验分布函数与假设的累积分布函数之间的最大正差异（D+），则为 +1；如果 KS 统计量是最大负差异（D-），则为 -1。

另请参见

ks_2samp，kstest

注意

有三种选项用于空假设和相应的备择假设，可以使用alternative参数进行选择。

• 双侧：零假设是两个分布相同，即 F(x)=G(x)对所有 x 成立；备择假设是它们不相同。

• 更少：零假设是对所有 x，F(x) >= G(x)成立；备择假设是对至少一个 x，F(x) < G(x)成立。

• 更大：零假设是对所有 x，F(x) <= G(x)成立；备择假设是对至少一个 x，F(x) > G(x)成立。

注意备择假设描述的是底层分布的CDF，而不是观察值。例如，假设 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.ks_1samp(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.ks_1samp(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.ks_1samp(x, stats.norm.cdf, alternative='less')
KstestResult(statistic=0.17482387821055168, pvalue=0.001913921057766743) 

而且，由于 p 值小于我们的阈值，我们拒绝零假设，支持备择假设。

scipy.stats.goodness_of_fit

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

scipy.stats.goodness_of_fit(dist, data, *, known_params=None, fit_params=None, guessed_params=None, statistic='ad', n_mc_samples=9999, random_state=None)

执行一个比较数据与分布族的拟合优度检验。

给定分布族和数据，执行空假设检验，即数据是否来自该族分布的检验。可以指定分布的任何已知参数。分布的剩余参数将拟合到数据中，并相应地计算检验的 p 值。提供了几种比较分布与数据的统计量。

参数：

distscipy.stats.rv_continuous

代表空假设下的分布族的对象。

data1D array_like

要测试的有限未经审查的数据。

known_paramsdict，可选

包含已知分布参数名称-值对的字典。蒙特卡罗样本从假设的空假设分布中随机抽取这些参数值。在每个蒙特卡罗样本中，只有剩余的空假设分布族的未知参数被拟合到样本中；已知参数保持不变。如果所有分布族参数都已知，则省略将分布族拟合到每个样本的步骤。

fit_paramsdict，可选

包含已经拟合到数据的分布参数名称-值对的字典，例如使用scipy.stats.fit或dist的fit方法。蒙特卡罗样本从假设的空假设分布中抽取，使用这些指定的参数值。然而，在这些蒙特卡罗样本上，空假设分布族的这些以及所有其他未知参数在计算统计量之前被拟合。

guessed_paramsdict，可选

包含已经猜测的分布参数名称-值对的字典。这些参数始终被视为自由参数，并且被拟合到提供的data以及从空假设分布中抽取的蒙特卡罗样本中。这些guessed_params的目的是作为数值拟合过程的初始值使用。

statistic{“ad”, “ks”, “cvm”, “filliben”}，可选

用于将数据与分布进行比较的统计量，在将分布族的未知参数拟合到数据之后进行。Anderson-Darling（“ad”）[1]、Kolmogorov-Smirnov（“ks”）[1]、Cramer-von Mises（“cvm”）[1] 和 Filliben（“filliben”）[7] 统计量可用。

n_mc_samplesint，默认值：9999

从零假设分布中抽取的蒙特卡洛样本数量。每个样本的样本量与给定的data相同。

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

numpy.random.RandomState，可选

用于生成蒙特卡洛样本的伪随机数生成器状态。

如果random_state为None（默认），则使用numpy.random.RandomState单例。如果random_state为整数，则使用新的RandomState实例，并使用random_state作为种子。如果random_state已经是Generator或RandomState实例，则使用提供的实例。

返回：

resGoodnessOfFitResult

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

fit_resultFitResult

一个表示提供的dist与data拟合情况的对象。此对象包括完全定义零假设分布的分布族参数值，即从中抽取蒙特卡洛样本的分布。

statisticfloat

比较提供的data与零假设分布的统计量值。

pvaluefloat

零分布中具有与提供的data的统计量值至少一样极端的元素的比例。

null_distributionndarray

每个从零假设分布抽取的蒙特卡洛样本的统计量值。

注意事项

这是一种广义的蒙特卡洛拟合优度检验过程，其特殊情况对应于各种 Anderson-Darling 测试、Lilliefors 测试等。该测试在文献[2]、[3]和[4]中被描述为参数化自举检验。这是一个蒙特卡洛检验，其中从数据中估计了用于抽取样本的分布的参数。我们在以下描述中使用“蒙特卡洛”而不是“参数化自举”，以避免与更熟悉的非参数化自举混淆，并描述了测试的执行方式。

传统的拟合优度检验

传统上，对应于固定的显著性水平集的临界值是使用蒙特卡洛方法预先计算的。用户通过仅计算他们观察到的数据的测试统计值并将此值与表格化的临界值进行比较来执行测试。这种做法不太灵活，因为并非所有分布和已知和未知参数值的组合都有可用的表格。当从有限的表格数据插值临界值以与用户的样本大小和拟合参数值对应时，结果可能不准确。为了克服这些缺点，此函数允许用户执行适应其特定数据的蒙特卡洛试验。

算法概述

简言之，此例程执行以下步骤：

1. 将未知参数拟合到给定的数据，从而形成“零假设”分布，并计算此数据和分布对的统计量。
3. 从这个零假设分布中抽取随机样本。
5. 将未知参数拟合到每个随机样本。
7. 计算每个样本与拟合到样本的分布之间的统计量。
9. 将来自（1）的与数据相对应的统计值与来自（4）的随机样本的统计值进行比较。p 值是具有大于或等于观察数据的统计值的样本比例。

更详细地说，步骤如下。

首先，使用最大似然估计将指定的分布家族dist的任何未知参数拟合到提供的数据中。（一个例外是具有未知位置和尺度的正态分布：我们使用偏差校正标准差 np.std(data, ddof=1) 作为尺度，如[1]中推荐的那样。）这些参数的值指定了分布家族的特定成员，称为“零假设”分布，即从中数据在零假设下进行抽样的分布。统计量，它比较数据与分布之间的关系，计算在数据和零假设分布之间的。

接下来，从零假设分布中抽取许多（具体为n_mc_samples）新样本，每个样本包含与数据相同数量的观测值。将分布家族dist的所有未知参数适应于每个重新采样，并计算每个样本与其相应拟合分布之间的统计量。这些统计量值形成蒙特卡洛零分布（不要与上面的“零假设”分布混淆）。

测试的 p 值是蒙特卡洛零分布中统计值的比例，这些统计值至少与所提供的数据的统计值一样极端。更确切地说，p 值由以下公式给出：

[p = \frac{b + 1} {m + 1}]

其中 (b) 是蒙特卡洛空分布中的统计值数量，这些值大于或等于为 data 计算的统计值，(m) 是蒙特卡洛空分布中元素的数量 (n_mc_samples)。将分子和分母各加 (1) 可以理解为将与 data 相对应的统计值包括在空分布中，但更正式的解释见文献 [5]。

限制

由于必须对分布族的未知参数逐个拟合蒙特卡洛样本中的每一个，对于某些分布族而言，该测试可能非常缓慢；在 SciPy 中的大多数分布，通过数值优化进行分布拟合。

反模式

由于这个原因，可能会诱使将分布的参数（由用户预先拟合到 data）视为 known_params，因此在每个蒙特卡洛样本中拟合分布的需求被排除。这本质上是如何执行原始的科尔莫哥洛夫-斯米尔诺夫检验的。尽管这样的检验可以提供反对零假设的证据，但这样的检验在某种意义上是保守的，因为小的 p 值倾向于（极大地）高估发生第一类错误的概率（即虽然接受了零假设，但它是真实的），而检验的功效较低（即，即使零假设是错误的，也不太可能拒绝零假设）。这是因为蒙特卡洛样本不太可能与零假设分布以及 data 一致，从而增加了空分布中记录的统计值，使得更多的统计值超过 data 的统计值，从而增加了 p 值。

参考文献

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

M. A. Stephens (1974). “适合性检验的经验分布函数统计量及其比较。” 《美国统计学会杂志》第 69 卷，第 730-737 页。

[2]

W. Stute, W. G. Manteiga, 和 M. P. Quindimil (1993). “基于自举法的拟合优度检验。” 《Metrika》40.1: 243-256。

[3]

C. Genest 和 B Rémillard (2008). “半参数模型中拟合优度检验的参数自举法的有效性。” 《法国数学与统计学院概率与统计学年刊》44.6。

[4]

I. Kojadinovic 和 J. Yan (2012). “基于加权自举的拟合优度检验：大样本参数自举的快速替代。” 《加拿大统计杂志》40.3: 480-500。

[5]

B. Phipson 和 G. K. Smyth (2010). “排列 P 值不应为零：在随机抽样排列时计算确切 P 值。” 《统计应用于遗传学和分子生物学》9.1。

[6]

H. W. Lilliefors (1967). “关于均值和方差未知的正态分布的科尔莫哥洛夫-斯米尔诺夫检验。” 《美国统计学会杂志》62.318: 399-402。

[7]

Filliben, James J.，“用于正态性的概率图相关系数检验。” Technometrics 17.1（1975）：111-117。

示例

一个广为人知的用于检验数据是否来自特定分布的零假设的测试是科尔莫哥罗夫-斯米尔诺夫（KS）检验，在 SciPy 中可通过scipy.stats.ks_1samp找到。假设我们希望测试以下数据：

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> x = stats.uniform.rvs(size=75, random_state=rng) 

从正态分布中抽样。要执行 KS 检验，将比较观察数据的经验分布函数与（理论上的）正态分布的累积分布函数。当然，为了做到这一点，必须完全指定零假设下的正态分布。通常首先通过将分布的loc和scale参数拟合到观察数据，然后执行测试来完成此操作。

>>> loc, scale = np.mean(x), np.std(x, ddof=1)
>>> cdf = stats.norm(loc, scale).cdf
>>> stats.ks_1samp(x, cdf)
KstestResult(statistic=0.1119257570456813, pvalue=0.2827756409939257) 

KS 测试的一个优点是可以精确和高效地计算 p 值 - 在零假设下获得测试统计量值的概率，该值与从观察数据中获得的值一样极端。goodness_of_fit只能近似这些结果。

>>> known_params = {'loc': loc, 'scale': scale}
>>> res = stats.goodness_of_fit(stats.norm, x, known_params=known_params,
...                             statistic='ks', random_state=rng)
>>> res.statistic, res.pvalue
(0.1119257570456813, 0.2788) 

检验统计量完全匹配，但 p 值是通过形成“蒙特卡洛零分布”来估计的，即通过从scipy.stats.norm中提供的参数显式抽取随机样本，并计算每个统计量。至少与res.statistic一样极端的这些统计值的比例近似于通过scipy.stats.ks_1samp计算的精确 p 值。

然而，在许多情况下，我们更愿意仅测试数据是否来自正态分布族的任何成员之一，而不是特别来自具有拟合到观察样本的位置和比例的正态分布。在这种情况下，Lilliefors [6]认为 KS 检验过于保守（即 p 值夸大了拒绝真空假设的实际概率），因此缺乏功效 - 即在真空假设实际为假时拒绝真空假设的能力。实际上，我们的 p 值约为 0.28，这远远大于在任何常见显著性水平下拒绝零假设的实际概率。

考虑为什么会这样。请注意，在上述 KS 检验中，统计量始终将数据与拟合到观察数据的正态分布的累积分布函数进行比较。这倾向于降低观察数据的统计量值，但在计算其他样本的统计量时（例如我们随机抽取的样本以形成蒙特卡罗零分布时），这种方式是“不公平”的。我们可以很容易地进行修正：每当我们计算样本的 KS 统计量时，我们使用拟合到该样本的正态分布的累积分布函数。在这种情况下，零分布未经精确计算，通常是使用上述的蒙特卡罗方法来近似的。这就是 goodness_of_fit 突出表现的地方。

>>> res = stats.goodness_of_fit(stats.norm, x, statistic='ks',
...                             random_state=rng)
>>> res.statistic, res.pvalue
(0.1119257570456813, 0.0196) 

实际上，这个 p 值要小得多，足够小以（正确地）在常见的显著水平下拒绝零假设，包括 5% 和 2.5%。

然而，KS 统计量对所有与正态分布偏差不是很敏感。KS 统计量最初的优势在于能够理论上计算零分布，但现在我们可以通过计算来近似零分布，可以使用更敏感的统计量 - 从而得到更高的检验力。Anderson-Darling 统计量 [1] 倾向于更为敏感，已经使用蒙特卡罗方法为各种显著水平和样本大小制表了此统计量的临界值。

>>> res = stats.anderson(x, 'norm')
>>> print(res.statistic)
1.2139573337497467
>>> print(res.critical_values)
[0.549 0.625 0.75  0.875 1.041]
>>> print(res.significance_level)
[15\.  10\.   5\.   2.5  1\. ] 

在这里，统计量的观察值超过了对应于 1% 显著水平的临界值。这告诉我们观察数据的 p 值小于 1%，但确切值是多少？我们可以从这些（已经插值的）值中插值，但 goodness_of_fit 可以直接估计它。

>>> res = stats.goodness_of_fit(stats.norm, x, statistic='ad',
...                             random_state=rng)
>>> res.statistic, res.pvalue
(1.2139573337497467, 0.0034) 

另一个优势是使用 goodness_of_fit 不受限于特定分布或已知参数与需从数据中估计参数的条件。相反，goodness_of_fit 可以对任何具有足够快速和可靠 fit 方法的分布相对快速地估计 p 值。例如，在这里我们使用 Cramer-von Mises 统计量对具有已知位置和未知尺度的 Rayleigh 分布进行拟合优度检验。

>>> rng = np.random.default_rng()
>>> x = stats.chi(df=2.2, loc=0, scale=2).rvs(size=1000, random_state=rng)
>>> res = stats.goodness_of_fit(stats.rayleigh, x, statistic='cvm',
...                             known_params={'loc': 0}, random_state=rng) 

这个过程执行起来非常快速，但是为了检查 fit 方法的可靠性，我们应该检查拟合结果。

>>> res.fit_result  # location is as specified, and scale is reasonable
 params: FitParams(loc=0.0, scale=2.1026719844231243)
 success: True
 message: 'The fit was performed successfully.'
>>> import matplotlib.pyplot as plt  # matplotlib must be installed to plot
>>> res.fit_result.plot()
>>> plt.show() 

如果分布未能尽可能地拟合观察数据，测试可能无法控制类型 I 错误率，即在零假设为真时拒绝零假设的概率。

我们还应该寻找零分布中的极端异常值，这些异常值可能是由于不可靠的拟合导致的。这些异常值不一定会使结果无效，但它们倾向于降低检验的功效。

>>> _, ax = plt.subplots()
>>> ax.hist(np.log10(res.null_distribution))
>>> ax.set_xlabel("log10 of CVM statistic under the null hypothesis")
>>> ax.set_ylabel("Frequency")
>>> ax.set_title("Histogram of the Monte Carlo null distribution")
>>> plt.show() 

这个图看起来令人放心。

如果 fit 方法可靠运行，并且测试统计量的分布对拟合参数的值不是特别敏感，那么由 goodness_of_fit 提供的 p 值预计会是一个很好的近似。

>>> res.statistic, res.pvalue
(0.2231991510248692, 0.0525) 

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_obs和f_exp在其上应用测试。如果 axis 为 None，则将f_obs中的所有值视为单个数据集。默认为 0。

返回：

res：Power_divergenceResult

包含属性的对象：

statisticfloat 或 ndarray

卡方检验统计量。如果 axis 为 None 或f_obs和f_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_obs 和 f_exp 也会进行广播。在下面的例子中，f_obs 的形状为 (6,)，f_exp 的形状为 (2, 6)，因此广播 f_obs 和 f_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_obs和f_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_obs和f_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_obs 或 f_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 检验](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_obs 和 f_exp 也在广播中使用。在下面的例子中，f_obs 的形状为 (6,)，f_exp 的形状为 (2, 6)，因此广播 f_obs 和 f_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’, ‘less’, ‘greater’}，可选

定义备择假设。有以下选项可用（默认为‘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的形式返回，包含low和high字段。

从 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 的元素在确切算术时被分配不同的排名，即使它们会因精确算术而绑定。与分开传递 x 和 y 不同，考虑计算差异 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 表示配对样本之间的差异：如果同时提供 x 和 y，则 d = x - y，否则 d = x。

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

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

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

method{“auto”, “exact”, “approx”}，可选

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

axis整数或 None，默认为 0

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

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

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

• propagate：如果轴切片（例如行）中存在 NaN，则输出的相应条目将为 NaN。

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

• raise：如果存在 NaN，则会引发 ValueError。

keepdimsbool，默认为 False

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

返回：

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

statistic类似数组

如果 alternative 是“双侧”，则是差异排名之和（无论是正还是负）。否则是正差异的排名之和。

pvalue类似数组

测试的 p 值取决于 alternative 和 method。

zstatistic类似数组

当 method = 'approx' 时，这是标准化的 z 统计量：

z = (T - mn - d) / se 

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

另请参阅

kruskal，mannwhitneyu

备注

在以下内容中，让 d 表示成对样本之间的差异：如果提供了 x 和 y，则 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 值，参考文献中并无明确共识。无论如何，这是当使用wilcoxon和pymethod='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接受等长的两个样本，计算配对元素之间的差异，然后进行测试。考虑样本 x 和 y：

>>> 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=None 且 x 是一个 2x2 数组的情况下，linregress(x) 等同于 linregress(x[0], x[1])。

alternative{‘two-sided’, ‘less’, ‘greater’}, 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，具有字段 slope、intercept、rvalue、pvalue 和 stderr，因此可以继续编写：

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() 

计算斜率和截距的 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’, ‘greater’, ‘less’}，可选

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

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

• ‘less’: 相关性为负（小于零）

• ‘greater’: 相关性为正（大于零）

自版本 1.9.0 新增。

methodResamplingMethod，可选

定义了计算 p 值的方法。如果 method 是 PermutationMethod/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的形式返回，字段为low和high。如果未提供method，则使用 Fisher 变换计算置信区间[1]。如果method是BootstrapMethod的一个实例，则使用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{‘propagate’, ‘raise’, ‘omit’}，可选

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

• ‘propagate’：返回 NaN

• ‘raise’：抛出错误

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

alternative{‘two-sided’, ‘less’, ‘greater’}，可选

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

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

• ‘less’：相关性为负（小于零）

• ‘greater’：相关性为正（大于零）

新版本为 1.7.0。

返回：

resSignificanceResult

一个包含属性的对象：

statisticfloat 或 ndarray（2-D 方阵）

Spearman 相关系数矩阵或相关系数（如果仅给出 2 个变量作为参数）。相关系数矩阵是方阵，其长度等于a和b组合后的总变量数（列或行）。

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]的数据，研究了不健康人类肝脏中游离脯氨酸（一种氨基酸）和总胶原（经常存在于结缔组织中的蛋白质）之间的关系。

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

>>> 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() 

比较通过 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() 

>>> 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() 

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

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

>>> 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() 

>>> 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{‘propagate’, ‘raise’, ‘omit’}，可选

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

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

method{‘auto’, ‘asymptotic’, ‘exact’}，可选

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

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

variant{‘b’, ‘c’}，可选

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

alternative{‘two-sided’, ‘less’, ‘greater’}，可选

定义备择假设。默认为‘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中的绑定数。如果同一对在x和y中都有绑定，则不会添加到 T 或 U 中。n 是样本的总数，m 是在x或y中较小的唯一值的数量。

参考文献

[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] 的以下数据，该研究了不健康的人类肝脏中自由脯氨酸（一种氨基酸）与总胶原（一种经常在结缔组织中找到的蛋白质）之间的关系。

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

>>> 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() 

比较通过 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() 

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

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

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

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

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

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

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

>>> 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() 

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

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