SciPy 1.12 中文文档（六十五）

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

scipy.stats.levene

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

scipy.stats.levene(*samples, center='median', proportiontocut=0.05, axis=0, nan_policy='propagate', keepdims=False)

执行列文氏检验以检验方差是否相等。

列文氏检验检验的是所有输入样本来自具有相等方差的总体的零假设。列文氏检验是巴特利特检验 bartlett 在存在明显偏离正态分布时的替代方法。

参数：

**sample1, sample2, …**array_like

样本数据，可能长度不同。只接受一维样本。

center{‘mean’, ‘median’, ‘trimmed’}, optional

在测试中使用数据的哪个函数。默认为 'median'。

proportiontocutfloat, optional

当 center 为 ‘trimmed’ 时，这给出了要从每端裁剪的数据点的比例。（见 scipy.stats.trim_mean.）默认为 0.05。

axisint or None, default: 0

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

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

定义如何处理输入的 NaN。

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

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

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

keepdimsbool, default: False

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

返回：

statisticfloat

检验统计量。

pvaluefloat

测试的 p 值。

另请参阅

fligner

正态样本中 k 个方差的非参数检验

bartlett

正态样本中 k 个方差的参数检验

注意

列文氏检验有三种变体。各种可能性及其建议的用法如下：

• ‘median’ : 建议用于偏斜（非正态）分布>
• ‘mean’ : 建议用于对称的，中等尾部的分布。
• ‘trimmed’ : 建议用于重尾分布。

Levene 的测试版本使用了均值，在 Levene 的原始文章中提出（[2]），而中位数和修剪均值由 Brown 和 Forsythe 研究（[3]），有时也称为 Brown-Forsythe 测试。

从 SciPy 1.9 开始，不推荐使用 np.matrix 输入，计算前会将其转换为 np.ndarray。在这种情况下，输出将是相应形状的标量或 np.ndarray，而不是二维 np.matrix。类似地，虽然掩码数组的掩码元素被忽略，输出将是标量或 np.ndarray，而不是具有 mask=False 的掩码数组。

参考文献

[1]

www.itl.nist.gov/div898/handbook/eda/section3/eda35a.htm

[2]

Levene, H.（1960）。在 Harold Hotelling 的《概率和统计的贡献：向》中，I. Olkin 等人编辑，斯坦福大学出版社，278-292 页。

[3]

Brown, M. B. 和 Forsythe, A. B.（1974），《美国统计协会杂志》，69，364-367

[4]

C.I. BLISS（1952），《生物测定统计学：特别参考维生素》，pp 499-503，DOI:10.1016/C2013-0-12584-6。

[5]

B. Phipson 和 G. K. Smyth。“当置换随机抽取时，置换 p 值不应为零：计算确切 p 值。”《统计应用于遗传学和分子生物学》，9.1（2010）。

[6]

Ludbrook, J. 和 Dudley, H.（1998）。《为什么在生物医学研究中，置换检验优于 t 检验和 F 检验》。《美国统计学家》，52（2），127-132。

例子

在 [4] 中，研究了维生素 C 对豚鼠牙齿生长的影响。在控制研究中，60 名受试者被分为小剂量、中剂量和大剂量组，分别每天服用 0.5、1.0 和 2.0 毫克的维生素 C。42 天后，测量了牙齿生长情况。

下面的 small_dose、medium_dose 和 large_dose 数组记录了三组牙齿生长的微米测量值。

>>> import numpy as np
>>> small_dose = np.array([
...     4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
...     15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7
... ])
>>> medium_dose = np.array([
...     16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
...     19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3
... ])
>>> large_dose = np.array([
...     23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
...     25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23
... ]) 

levene 统计量对样本间方差差异敏感。

>>> from scipy import stats
>>> res = stats.levene(small_dose, medium_dose, large_dose)
>>> res.statistic
0.6457341109631506 

当样本方差差异较大时，统计量的值往往较高。

我们可以通过比较统计量的观察值与零分布来测试组间方差的不等性：即在零假设下，三组总体方差相等的假设下得到的统计值分布。

对于这个测试，零分布如下图所示，遵循 F 分布。

>>> import matplotlib.pyplot as plt
>>> k, n = 3, 60   # number of samples, total number of observations
>>> dist = stats.f(dfn=k-1, dfd=n-k)
>>> val = np.linspace(0, 5, 100)
>>> pdf = dist.pdf(val)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def plot(ax):  # we'll reuse this
...     ax.plot(val, pdf, color='C0')
...     ax.set_title("Levene Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
...     ax.set_xlim(0, 5)
...     ax.set_ylim(0, 1)
>>> plot(ax)
>>> plt.show() 

比较由 p 值量化：即零分布中大于或等于观察统计值的比例。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> pvalue = dist.sf(res.statistic)
>>> annotation = (f'p-value={pvalue:.3f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (1.5, 0.22), (2.25, 0.3), arrowprops=props)
>>> i = val >= res.statistic
>>> ax.fill_between(val[i], y1=0, y2=pdf[i], color='C0')
>>> plt.show() 

>>> res.pvalue
0.5280694573759905 

如果 p 值很“小” - 也就是说，从具有相同方差的分布中抽样数据产生了如此极端的统计值的概率很低 - 这可能被视为反对零假设的证据，有利于替代假设：组的方差不相等。注意：

• 反之则不成立；也就是说，这个测试不能用来证明零假设。

• 被认为“小”的值的阈值应在分析数据之前做出选择，考虑到假阳性（错误拒绝零假设）和假阴性（未能拒绝虚假零假设）的风险[5]。

• 小的 p 值并不是对大效应的证据；相反，它们只能提供对“显著”效应的证据，这意味着在零假设下这些结果发生的可能性很小。

注意 F 分布提供了零分布的渐近近似。对于小样本，执行置换检验可能更合适：在零假设下，所有三个样本都是从同一总体中抽取的，每个测量值等可能地出现在三个样本中的任何一个。因此，我们可以通过在观察值随机分区的许多随机生成中计算统计量来形成随机化的零分布。

>>> def statistic(*samples):
...     return stats.levene(*samples).statistic
>>> ref = stats.permutation_test(
...     (small_dose, medium_dose, large_dose), statistic,
...     permutation_type='independent', alternative='greater'
... )
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> bins = np.linspace(0, 5, 25)
>>> ax.hist(
...     ref.null_distribution, bins=bins, density=True, facecolor="C1"
... )
>>> ax.legend(['aymptotic approximation\n(many observations)',
...            'randomized null distribution'])
>>> plot(ax)
>>> plt.show() 

>>> ref.pvalue  # randomized test p-value
0.4559  # may vary 

注意，这里计算的 p 值与levene上返回的渐近近似之间存在显著分歧。可以从置换检验中严格推导出的统计推断有限；尽管如此，在许多情况下，这可能是首选方法[6]。

下面是另一个一般性示例，其中零假设将被拒绝。

测试列表a、b和c是否来自具有相等方差的总体。

>>> a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
>>> b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
>>> c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
>>> stat, p = stats.levene(a, b, c)
>>> p
0.002431505967249681 

小的 p 值表明这些总体的方差不相等。

鉴于b的样本方差远大于a和c的样本方差，这并不令人惊讶：

>>> [np.var(x, ddof=1) for x in [a, b, c]]
[0.007054444444444413, 0.13073888888888888, 0.008890000000000002] 

scipy.stats.bartlett

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

scipy.stats.bartlett(*samples, axis=0, nan_policy='propagate', keepdims=False)

执行 Bartlett 的等方差性检验。

Bartlett 测试检验所有输入样本是否来自具有相等方差的总体的零假设。对于来自显著非正态分布的样本，Levene 的测试 levene 更为稳健。

参数：

**sample1, sample2, …**array_like

样本数据的数组。仅接受 1d 数组，它们可能有不同的长度。

轴int 或 None，默认为 0

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

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

定义如何处理输入 NaN。

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

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

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

keepdimsbool，默认为 False

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

返回：

统计量float

检验统计量。

pvaluefloat

检验的 p 值。

另请参阅

fligner

一种用于 k 个方差相等的非参数检验

levene

一种用于 k 个方差相等的稳健参数检验

笔记

Conover 等人（1981）通过大量模拟研究了许多现有的参数和非参数检验，并得出结论，Fligner 和 Killeen（1976）以及 Levene（1960）提出的检验在正态性偏差和功效方面似乎更为优越（[3]）。

自 SciPy 1.9 开始，将 np.matrix 输入（不建议新代码使用）转换为 np.ndarray 后再执行计算。在这种情况下，输出将是适当形状的标量或 np.ndarray，而不是 2D np.matrix。同样，虽然忽略掩码数组的屏蔽元素，但输出将是标量或 np.ndarray，而不是具有 mask=False 的屏蔽数组。

参考文献

[1]

www.itl.nist.gov/div898/handbook/eda/section3/eda357.htm

[2]

Snedecor, George W. and Cochran, William G. (1989), 统计方法，第八版，爱荷华州立大学出版社。

[3]

Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and Hypothesis Testing based on Quadratic Inference Function. Technical Report #99-03, Center for Likelihood Studies, Pennsylvania State University.

[4]

Bartlett, M. S. (1937). Sufficiency and Statistical Tests 的特性。 伦敦皇家学会会议录 A 系列，数学和物理科学，第 160 卷，第 901 号，第 268-282 页。

[5]

C.I. BLISS (1952), 生物测定统计学：特别参考维生素，第 499-503 页，DOI:10.1016/C2013-0-12584-6。

[6]

B. Phipson and G. K. Smyth. “置换 P 值永远不应为零：当置换随机抽取时计算确切的 P 值。” 遗传学和分子生物学中的统计应用 9.1（2010 年）。

[7]

Ludbrook, J., & Dudley, H. (1998). 为什么在生物医学研究中，置换检验优于 t 和 F 检验。 美国统计学家，52（2），127-132。

例子

在[5]中，研究了维生素 C 对豚鼠牙齿生长的影响。 在一项对照研究中，60 名受试者分为小剂量，中剂量和大剂量组，分别每天摄取 0.5、1.0 和 2.0 毫克维生素 C。 42 天后，测量了牙齿生长情况。

下面的small_dose，medium_dose和large_dose数组记录了三组的牙齿生长测量结果，单位是微米。

>>> import numpy as np
>>> small_dose = np.array([
...     4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
...     15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7
... ])
>>> medium_dose = np.array([
...     16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
...     19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3
... ])
>>> large_dose = np.array([
...     23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
...     25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23
... ]) 

bartlett 统计量对样本之间的方差差异敏感。

>>> from scipy import stats
>>> res = stats.bartlett(small_dose, medium_dose, large_dose)
>>> res.statistic
0.6654670663030519 

当方差差异很大时，统计量的值往往较高。

我们可以通过比较统计量观察值与零假设下的零分布来检验组间方差的不平等：假设三组的总体方差相等的零假设下的统计量值的分布。

对于此测试，零分布遵循如下卡方分布。

>>> import matplotlib.pyplot as plt
>>> k = 3  # number of samples
>>> dist = stats.chi2(df=k-1)
>>> val = np.linspace(0, 5, 100)
>>> pdf = dist.pdf(val)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def plot(ax):  # we'll reuse this
...     ax.plot(val, pdf, color='C0')
...     ax.set_title("Bartlett Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
...     ax.set_xlim(0, 5)
...     ax.set_ylim(0, 1)
>>> plot(ax)
>>> plt.show() 

比较通过 p 值来量化：在零分布中大于或等于统计量观察值的值的比例。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> pvalue = dist.sf(res.statistic)
>>> annotation = (f'p-value={pvalue:.3f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (1.5, 0.22), (2.25, 0.3), arrowprops=props)
>>> i = val >= res.statistic
>>> ax.fill_between(val[i], y1=0, y2=pdf[i], color='C0')
>>> plt.show() 

>>> res.pvalue
0.71696121509966 

如果 p 值“很小” - 也就是说，从具有相同方差的分布中抽取数据产生这样极端统计值的概率很低 - 这可能被视为支持备择假设的证据：组的方差不相等。请注意：

• 逆不成立；也就是说，此测试不能用来支持零假设。

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

• 较小的 p 值并不是大效应的证据；相反，它们只能为“显著”效应提供证据，意味着在零假设下不太可能发生。

请注意，当观测值服从正态分布时，卡方分布提供零分布。对于从非正态总体中抽取的小样本，执行置换检验可能更为合适：在零假设下，所有三个样本均从同一总体中抽取，每个测量值在三个样本中被观察到的概率相等。因此，我们可以通过在观测值随机分割成三个样本的许多随机生成的分割中计算统计量来形成随机化的零分布。

>>> def statistic(*samples):
...     return stats.bartlett(*samples).statistic
>>> ref = stats.permutation_test(
...     (small_dose, medium_dose, large_dose), statistic,
...     permutation_type='independent', alternative='greater'
... )
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> bins = np.linspace(0, 5, 25)
>>> ax.hist(
...     ref.null_distribution, bins=bins, density=True, facecolor="C1"
... )
>>> ax.legend(['aymptotic approximation\n(many observations)',
...            'randomized null distribution'])
>>> plot(ax)
>>> plt.show() 

>>> ref.pvalue  # randomized test p-value
0.5387  # may vary 

请注意，此处计算的 p 值与bartlett提供的渐近近似存在显著差异。从置换检验中可以严格推导出的统计推断有限；尽管如此，在许多情况下，它们可能是首选方法[7]。

以下是另一个通用示例，拒绝零假设的情况。

检验列表a、b和c是否来自具有相等方差的总体。

>>> a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
>>> b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
>>> c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
>>> stat, p = stats.bartlett(a, b, c)
>>> p
1.1254782518834628e-05 

非常小的 p 值表明这些总体的方差不相等。

这并不令人意外，因为b的样本方差远大于a和c的样本方差：

>>> [np.var(x, ddof=1) for x in [a, b, c]]
[0.007054444444444413, 0.13073888888888888, 0.008890000000000002] 

scipy.stats.median_test

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

scipy.stats.median_test(*samples, ties='below', correction=True, lambda_=1, nan_policy='propagate')

执行 Mood 中位数检验。

检验两个或多个样本是否来自具有相同中位数的总体。

设 n = len(samples) 表示样本数。计算所有数据的“总中位数”，并通过将每个样本中的值分类为高于或低于总中位数来形成列联表。列联表与 correction 和 lambda_ 一起传递给 scipy.stats.chi2_contingency 计算检验统计量和 p 值。

参数：

**sample1, sample2, …**array_like

样本集。必须至少有两个样本。每个样本必须是包含至少一个值的一维序列。不要求样本具有相同的长度。

tiesstr，可选

确定在列联表中如何分类等于总中位数的值。该字符串必须是以下之一：

"below":
    Values equal to the grand median are counted as "below".
"above":
    Values equal to the grand median are counted as "above".
"ignore":
    Values equal to the grand median are not counted. 

默认为“below”。

correctionbool，可选

如果为 True，并且只有两个样本，则在计算与列联表相关的检验统计量时应用 Yates 修正。默认值为 True。

lambda_float 或 str，可选

默认情况下，在此检验中计算的统计量是 Pearson 卡方统计量。lambda_ 允许使用 Cressie-Read 功率差异族中的统计量。有关详细信息，请参阅 power_divergence。默认值为 1（Pearson 卡方统计量）。

nan_policy{‘propagate’, ‘raise’, ‘omit’}，可选

定义如何处理输入包含 NaN 时的情况。‘propagate’ 返回 NaN，‘raise’ 抛出错误，‘omit’ 在执行计算时忽略 NaN 值。默认为 ‘propagate’。

返回：

resMedianTestResult

包含属性的对象：

statisticfloat

检验统计量。返回的统计量由 lambda_ 决定。默认为 Pearson 卡方统计量。

pvaluefloat

检验的 p 值。

medianfloat

总中位数。

tablendarray

离散表。表的形状为 (2, n)，其中 n 是样本数。第一行保存大于总体中位数的值的计数，第二行保存小于总体中位数的值的计数。该表允许进一步分析，例如使用 scipy.stats.chi2_contingency 或者如果有两个样本，则使用 scipy.stats.fisher_exact 而无需重新计算表。如果 nan_policy 是 “propagate” 并且输入中存在 NaN，则 table 的返回值为 None。

请参阅

kruskal

对独立样本计算 Kruskal-Wallis H 检验。

mannwhitneyu

计算样本 x 和 y 的 Mann-Whitney 等级检验。

注释

新版本为 0.15.0。

参考文献

[1]

Mood, A. M., 统计理论导论。McGraw-Hill (1950), 第 394-399 页。

[2]

Zar, J. H., 生物统计分析, 第 5 版。Prentice Hall (2010). 见第 8.12 和 10.15 节。

示例

一个生物学家进行了一项实验，其中有三组植物。第 1 组有 16 棵植物，第 2 组有 15 棵植物，第 3 组有 17 棵植物。每棵植物产生若干颗种子。每组的种子计数如下：

Group 1: 10 14 14 18 20 22 24 25 31 31 32 39 43 43 48 49
Group 2: 28 30 31 33 34 35 36 40 44 55 57 61 91 92 99
Group 3:  0  3  9 22 23 25 25 33 34 34 40 45 46 48 62 67 84 

下面的代码将 Mood 的中位数检验应用于这些样本。

>>> g1 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
>>> g2 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]
>>> g3 = [0, 3, 9, 22, 23, 25, 25, 33, 34, 34, 40, 45, 46, 48, 62, 67, 84]
>>> from scipy.stats import median_test
>>> res = median_test(g1, g2, g3) 

中位数是

>>> res.median
34.0 

并且离散表是

>>> res.table
array([[ 5, 10,  7],
 [11,  5, 10]]) 

p 太大，无法得出中位数不相同的结论：

>>> res.pvalue
0.12609082774093244 

“G 检验”可以通过将 lambda_="log-likelihood" 传递给 median_test 来执行。

>>> res = median_test(g1, g2, g3, lambda_="log-likelihood")
>>> res.pvalue
0.12224779737117837 

中位数在数据中出现多次，例如，如果使用 ties="above"，则会得到不同的结果：

>>> res = median_test(g1, g2, g3, ties="above")
>>> res.pvalue
0.063873276069553273 

>>> res.table
array([[ 5, 11,  9],
 [11,  4,  8]]) 

此示例说明，如果数据集不大并且存在与中位数相等的值，则 p 值可能对 ties 的选择敏感。

scipy.stats.friedmanchisquare

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

scipy.stats.friedmanchisquare(*samples)

计算重复样本的 Friedman 检验。

Friedman 检验用于检验同一群体的重复样本具有相同分布的零假设。通常用于检验以不同方式获得的样本的一致性。例如，如果在同一组个体上使用了两种采样技术，可以使用 Friedman 检验来确定这两种采样技术是否一致。

参数：

**sample1, sample2, sample3…**array_like

观察数组。所有数组的元素数量必须相同。至少需要提供三个样本。

返回：

statisticfloat

测试统计量，考虑并校正并列数。

pvaluefloat

假设测试统计量服从卡方分布时的相关 p 值。

注释

由于假设测试统计量服从卡方分布，所以 p 值仅在 n > 10 且重复样本超过 6 次时才可靠。

参考文献

[1]

en.wikipedia.org/wiki/Friedman_test

[2]

P. Sprent and N.C. Smeeton，《应用非参数统计方法，第三版》。第六章，第 6.3.2 节。

示例

在[2]中，对一组七名学生进行了运动前、运动后立即以及运动后 5 分钟的脉搏率（每分钟）。是否有证据表明这三个场合的脉搏率相似？

我们首先提出零假设 (H_0)：

这三个场合的脉搏率相同。

让我们用 Friedman 检验来评估这一假设的合理性。

>>> from scipy.stats import friedmanchisquare
>>> before = [72, 96, 88, 92, 74, 76, 82]
>>> immediately_after = [120, 120, 132, 120, 101, 96, 112]
>>> five_min_after = [76, 95, 104, 96, 84, 72, 76]
>>> res = friedmanchisquare(before, immediately_after, five_min_after)
>>> res.statistic
10.57142857142857
>>> res.pvalue
0.005063414171757498 

使用 5%的显著性水平，我们会拒绝零假设，支持备择假设：“这三个场合的脉搏率不同”。

scipy.stats.anderson_ksamp

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

scipy.stats.anderson_ksamp(samples, midrank=True, *, method=None)

k-样本 Anderson-Darling 测试。

k-样本 Anderson-Darling 测试是单样本 Anderson-Darling 测试的修改。它测试零假设，即 k 个样本来自相同的总体，而无需指定该总体的分布函数。临界值取决于样本数量。

参数：

samples1-D 数组序列

数据样本数组中的数组。

midrankbool, 可选

Anderson-Darling 测试的类型，计算的是。默认情况下（True），是适用于连续和离散总体的中位秩测试。如果为 False，则使用右侧经验分布。

methodPermutationMethod, 可选

定义用于计算 p 值的方法。如果 method 是 PermutationMethod 的一个实例，则使用提供的配置选项和其他适当的设置计算 p 值。否则，p 值从表格化的值中插值。

返回：

resAnderson_ksampResult

一个包含属性的对象：

statisticfloat

规范化 k-样本 Anderson-Darling 测试统计量。

critical_valuesarray

显著水平为 25%，10%，5%，2.5%，1%，0.5%，0.1% 的临界值。

pvaluefloat

测试的近似 p 值。如果未提供 method，该值被截断 / 上限为 0.1% / 25%。

引发：

ValueError

如果提供的样本少于 2 个，一个样本为空，或样本中没有不同的观测值。

另请参见

ks_2samp

2 样本 Kolmogorov-Smirnov 测试

anderson

1 样本 Anderson-Darling 测试

注意

[1] 定义了 k-样本 Anderson-Darling 测试的三个版本：一个用于连续分布，两个用于可能发生样本之间的绑定的离散分布，在这些版本中默认情况下使用中位秩经验分布函数。此例程的默认值是计算基于中位秩经验分布函数的版本。此测试适用于连续和离散数据。如果将 midrank 设置为 False，则用于离散数据的右侧经验分布。

对应于显著性水平从 0.01 到 0.25 的临界值来自[1]。p 值被限制在 0.1% / 25%之间。由于未来版本可能扩展临界值的范围，建议不要测试p == 0.25，而是测试p >= 0.25（下限类似处理）。

新功能在版本 0.14.0 中引入。

参考文献

[1] (1,2,3)

Scholz, F. W 和 Stephens, M. A.（1987），K-Sample Anderson-Darling Tests，美国统计协会杂志，第 82 卷，第 918-924 页。

示例

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> res = stats.anderson_ksamp([rng.normal(size=50),
... rng.normal(loc=0.5, size=30)])
>>> res.statistic, res.pvalue
(1.974403288713695, 0.04991293614572478)
>>> res.critical_values
array([0.325, 1.226, 1.961, 2.718, 3.752, 4.592, 6.546]) 

由于返回的检验值大于 5%的临界值（1.961），可以在 5%水平上拒绝两个随机样本来自同一分布的零假设，但在 2.5%水平上不能。插值给出了约为 4.99%的近似 p 值。

>>> samples = [rng.normal(size=50), rng.normal(size=30),
...            rng.normal(size=20)]
>>> res = stats.anderson_ksamp(samples)
>>> res.statistic, res.pvalue
(-0.29103725200789504, 0.25)
>>> res.critical_values
array([ 0.44925884,  1.3052767 ,  1.9434184 ,  2.57696569,  3.41634856,
 4.07210043, 5.56419101]) 

对于来自相同分布的三个样本，无法拒绝零假设。报告的 p 值（25%）已被限制，可能不太准确（因为它对应于值 0.449，而统计量为-0.291）。

在 p 值被限制或样本量较小时，置换检验可能更准确。

>>> method = stats.PermutationMethod(n_resamples=9999, random_state=rng)
>>> res = stats.anderson_ksamp(samples, method=method)
>>> res.pvalue
0.5254 

scipy.stats.monte_carlo_test

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

scipy.stats.monte_carlo_test(data, rvs, statistic, *, vectorized=None, n_resamples=9999, batch=None, alternative='two-sided', axis=0)

执行蒙特卡洛假设检验。

data包含一个样本或一个或多个样本的序列。rvs指定了在空假设下data中样本的分布。给定data的statistic的值与使用rvs生成的n_resamples组样本集的统计量值进行比较。这给出了 p 值，即在空假设下观察到测试统计量的如此极端值的概率。

参数：

data数组或者数组序列

一组或者一系列观测值的数组。

rvs可调用函数或者可调用函数的元组

在空假设下生成随机变量的一个可调用函数或者一个可调用函数的序列。rvs的每个元素必须是一个接受关键字参数size的可调用函数（例如rvs(size=(m, n))），并返回该形状的 N 维数组样本。如果rvs是一个序列，则rvs中的可调用函数的数量必须与data中的样本数量匹配，即len(rvs) == len(data)。如果rvs是一个单一的可调用函数，则data被视为单个样本。

statistic可调用函数

要计算假设检验的 p 值的统计量。statistic必须是一个可调用函数，接受一个样本（例如statistic(sample)）或者len(rvs)个单独的样本（例如如果rvs包含两个可调用函数且data包含两个样本，则为statistic(samples1, sample2)），并返回相应的统计量。如果设置了vectorized为True，statistic还必须接受一个关键字参数axis，并且被向量化以计算沿着data中提供的axis的样本的统计量。

vectorized布尔值，可选

如果设置vectorized为False，则statistic不会传递关键字参数axis，并且预期仅计算 1D 样本的统计量。如果为True，则在传递 ND 样本数组时，statistic将传递关键字参数axis并且预期沿着axis计算统计量。如果为None（默认），如果statistic有参数axis，则vectorized将设置为True。使用向量化统计量通常可以减少计算时间。

n_resamples整数，默认值：9999

从rvs的每个可调用函数中抽取的样本数量。等效地，作为蒙特卡洛空假设使用的统计值的数量。

batch整数，可选

每次对statistic的调用中要处理的蒙特卡洛样本数。内存使用为 O(batch * sample.size[axis] )。默认为None，此时batch等于n_resamples。

alternative{‘双侧’, ‘小于’, ‘大于’}

要计算 p 值的备择假设。对于每个备择假设，p 值定义如下。

• 'greater' : 空分布中大于或等于观察到的检验统计量值的百分比。

• 'less' : 空分布中小于或等于观察到的检验统计量值的百分比。

• 'two-sided' : 上述 p 值中较小者的两倍。

axisint，默认值：0

data的轴（或data中的每个样本）用于计算统计量。

返回：

resMonteCarloTestResult

一个带有属性的对象：

statisticfloat or ndarray

观察数据的检验统计量。

pvaluefloat or ndarray

给定备选假设的 p 值。

null_distributionndarray

测试统计量在零假设下生成的值。

参考文献

[1]

B. Phipson 和 G. K. Smyth. “Permutation P-values Should Never Be Zero: Calculating Exact P-values When Permutations Are Randomly Drawn.” 统计遗传学和分子生物学应用 9.1 (2010).

示例

假设我们希望检验一个小样本是否来自正态分布。我们决定使用样本的偏度作为检验统计量，并且我们将认为 p 值为 0.05 是统计学显著的。

>>> import numpy as np
>>> from scipy import stats
>>> def statistic(x, axis):
...     return stats.skew(x, axis) 

收集数据后，我们计算检验统计量的观察值。

>>> rng = np.random.default_rng()
>>> x = stats.skewnorm.rvs(a=1, size=50, random_state=rng)
>>> statistic(x, axis=0)
0.12457412450240658 

要确定如果样本是从正态分布中抽取的，则观察到偏度的极端值的观察概率，我们可以执行蒙特卡罗假设检验。该测试将从它们的正态分布中随机抽取许多样本，计算每个样本的偏度，并将我们原始的偏度与此分布进行比较，以确定一个近似的 p 值。

>>> from scipy.stats import monte_carlo_test
>>> # because our statistic is vectorized, we pass `vectorized=True`
>>> rvs = lambda size: stats.norm.rvs(size=size, random_state=rng)
>>> res = monte_carlo_test(x, rvs, statistic, vectorized=True)
>>> print(res.statistic)
0.12457412450240658
>>> print(res.pvalue)
0.7012 

在零假设下，获得一个小于或等于观察值的检验统计量的概率约为 70%。这比我们选择的 5%阈值要大，因此我们不能将其视为反对零假设的显著证据。

请注意，这个 p 值基本上与scipy.stats.skewtest的 p 值相匹配，后者依赖于基于样本偏度的渐近分布的检验统计量。

>>> stats.skewtest(x).pvalue
0.6892046027110614 

对于小样本量，这个渐近逼近是无效的，但可以使用monte_carlo_test来处理任何大小的样本。

>>> x = stats.skewnorm.rvs(a=1, size=7, random_state=rng)
>>> # stats.skewtest(x) would produce an error due to small sample
>>> res = monte_carlo_test(x, rvs, statistic, vectorized=True) 

提供测试统计量的蒙特卡罗分布以便进一步研究。

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> ax.hist(res.null_distribution, bins=50)
>>> ax.set_title("Monte Carlo distribution of test statistic")
>>> ax.set_xlabel("Value of Statistic")
>>> ax.set_ylabel("Frequency")
>>> plt.show() 

scipy.stats.permutation_test

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

scipy.stats.permutation_test(data, statistic, *, permutation_type='independent', vectorized=None, n_resamples=9999, batch=None, alternative='two-sided', axis=0, random_state=None)

在提供的数据上对给定统计量进行置换检验。

对于独立样本统计量，零假设是数据是从相同分布中随机抽取的。对于配对样本统计量，可以测试两个零假设：数据被随机配对，或者数据被随机分配到样本中。

参数：

data类数组的可迭代对象

包含样本的数组，每个样本都是一组观测值。样本数组的维度必须与广播兼容，除了 axis 外。

statistic可调用对象

用于计算假设检验的 p 值的统计量。statistic 必须是一个可调用的函数，接受样本作为单独的参数（例如 statistic(*data)），并返回结果统计量。如果设置了 vectorized 为 True，则 statistic 还必须接受一个关键字参数 axis 并进行向量化以沿着样本数组的提供的 axis 计算统计量。

permutation_type{'independent', 'samples', 'pairings'}，可选

要执行的置换类型，符合零假设的要求。前两种置换类型适用于配对样本统计量，其中所有样本包含相同数量的观测值，并且沿着 axis 具有相应索引的观测值被认为是配对的；第三种适用于独立样本统计量。

• 'samples'：观测值被分配到不同的样本，但与其他样本中相同的观测值保持配对。这种置换类型适用于配对样本假设检验，如威尔科克森符号秩检验和配对 t 检验。

• 'pairings'：观测值与不同的观测值配对，但它们仍然在同一样本内。这种置换类型适用于具有统计量如斯皮尔曼相关系数 (\rho)、肯德尔 (\tau) 和皮尔逊 (r) 的关联/相关性检验。

• 'independent'（默认）：观测值被分配到不同的样本中。样本可以包含不同数量的观测值。这种置换类型适用于独立样本假设检验，如曼-惠特尼 U 检验和独立样本 t 检验。

  请参阅下面的注释部分以获取有关置换类型更详细的描述。

vectorized布尔值，可选

如果将 vectorized 设置为 False，则不会传递关键字参数 axis 给 statistic，并且期望它仅为 1D 样本计算统计量。如果为 True，则在传递 ND 样本数组时，将传递关键字参数 axis 给 statistic 并且期望沿着 axis 计算统计量。如果为 None（默认），如果 axis 是 statistic 的参数，则 vectorized 将设置为 True。使用矢量化统计量通常可以减少计算时间。

n_resamplesint 或 np.inf，默认值：9999

用于近似空值分布的随机排列（重新取样）的数量。如果大于或等于不同排列的数量，则将计算精确的空值分布。注意，随着样本大小的增长，不同排列的数量会非常迅速地增加，因此仅对非常小的数据集适用精确测试。

batchint，可选

每次调用statistic时处理的排列数量。内存使用量为 O(batch * n )，其中 n 是所有样本的总大小，不管 vectorized 的值如何。默认为 None，此时 batch 是排列的数量。

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

用于计算 p 值的备择假设。对于每个备择假设，p 值的定义如下。

• 'greater'：空值分布中大于或等于测试统计量观察值的百分比。

• 'less'：空值分布中小于或等于测试统计量观察值的百分比。

• 'two-sided'（默认）：上述 p 值之一的两倍较小的值。

注意，随机化测试的 p 值是根据[2]和[3]中建议的保守（过估计）近似计算的，而不是建议的无偏估计器[4]。也就是说，在计算随机化空值分布中与测试统计量观察值一样极端的比例时，分子和分母的值都增加了一。这种调整的解释是，测试统计量的观察值总是作为随机化空值分布的一个元素。用于双边 p 值的约定不是普遍适用的；如果喜欢不同的定义，则返回观察到的测试统计量和空值分布。

axisint，默认值：0

(广播)样本的轴，用于计算统计量。如果样本具有不同维数，则在考虑 axis 之前，对具有较少维度的样本前置单例维度。

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 实例，则使用该实例。

返回：

resPermutationTestResult

具有以下属性的对象：

statisticfloat 或 ndarray

数据的观察检验统计量。

pvaluefloat 或 ndarray

给定备择假设的 p 值。

null_distributionndarray

在零假设下生成的检验统计量值。

注意事项

此函数支持的三种排列检验类型如下所述。

非配对统计量 (permutation_type='independent'):

与此排列类型相关联的零假设是，所有观察值都从相同的基础分布中抽取，并且它们被随机分配到一个样本中。

假设 data 包含两个样本；例如 a, b = data。当 1 < n_resamples < binom(n, k) 时，其中

• k 是 a 中观测值的数量，

• n 是 a 和 b 中观测值的总数，以及

• binom(n, k) 是二项式系数 (n 选择 k)，

数据被合并（串联），随机分配到第一或第二个样本，并计算统计量。此过程重复执行 permutation 次，生成零假设下统计量的分布。将原始数据的统计量与该分布进行比较，以确定 p 值。

当 n_resamples >= binom(n, k) 时，执行精确检验：数据在每种不同的方式下精确地一次性分配到样本中，并形成精确的零假设分布。请注意，对于给定数据在样本之间的分区方式，仅考虑数据在每个样本内的一种排序/排列。对于不依赖于数据顺序在样本内的统计量来说，这显著降低了计算成本，而不会影响零分布的形状（因为每个值的频率/计数受相同因素影响）。

对于 a = [a1, a2, a3, a4] 和 b = [b1, b2, b3]，此排列类型的示例是 x = [b3, a1, a2, b2] 和 y = [a4, b1, a3]。因为精确检验仅考虑数据在每个样本内的一种排序/排列，所以像 x = [b3, a1, b2, a2] 和 y = [a4, a3, b1] 这样的重新采样不被视为与上述示例不同。

permutation_type='independent' 不支持单样本统计量，但可应用于具有超过两个样本的统计量。在这种情况下，如果 n 是每个样本中观测值数量的数组，则不同分区的数量是：

np.prod([binom(sum(n[i:]), sum(n[i+1:])) for i in range(len(n)-1)]) 

配对统计量，排列配对 (permutation_type='pairings'):

与此置换类型相关的零假设是，每个样本内的观测值都来自相同的基础分布，并且与其他样本元素的配对是随机的。

假设 data 只包含一个样本；例如 a, = data，我们希望考虑将 a 的元素与第二个样本 b 的元素的所有可能配对。设 n 是 a 中的观测数，也必须等于 b 中的观测数。

当 1 < n_resamples < factorial(n) 时，对 a 中的元素进行随机置换。用户提供的统计量接受一个数据参数，例如 a_perm，并计算考虑 a_perm 和 b 的统计量。重复执行这一过程，permutation 次，生成零假设下统计量的分布。将原始数据的统计量与该分布进行比较，以确定 p 值。

当 n_resamples >= factorial(n) 时，执行精确检验：对 a 按每种不同方式精确置换一次。因此，对 a 和 b 之间的每个唯一配对样本计算统计量一次。

对于 a = [a1, a2, a3] 和 b = [b1, b2, b3]，这种置换类型的示例是 a_perm = [a3, a1, a2]，而 b 保持原始顺序。

permutation_type='pairings' 支持包含任意数量样本的 data，每个样本必须包含相同数量的观测值。data 中提供的所有样本都独立进行置换。因此，如果 m 是样本数，n 是每个样本中的观测数，则精确检验的置换数为：

factorial(n)**m 

请注意，如果例如双样本统计量并不直接依赖于提供观测值的顺序 - 只依赖于观测值的配对，那么在 data 中只需提供其中一个样本。这大大降低了计算成本，但不影响零分布的形状（因为每个值的频率/计数受相同因素影响）。

配对统计，样本置换 (permutation_type='samples'):

与此置换类型相关的零假设是，每对观测值都来自相同的基础分布，并且它们被分配到的样本是随机的。

假设 data 包含两个样本；例如 a, b = data。设 n 是 a 中的观测数，也必须等于 b 中的观测数。

当 1 < n_resamples < 2**n 时，对 a 和 b 中的元素进行随机交换（保持它们的配对关系），并计算统计量。重复执行这一过程，permutation 次，生成零假设下统计量的分布。将原始数据的统计量与该分布进行比较，以确定 p 值。

当 n_resamples >= 2**n 时，执行精确检验：观察值被准确地分配到两个样本中的每一种不同方式（同时保持配对）一次。

对于 a = [a1, a2, a3] 和 b = [b1, b2, b3]，这种排列类型的一个示例是 x = [b1, a2, b3] 和 y = [a1, b2, a3]。

permutation_type='samples' 支持 data 包含任意数量的样本，每个样本必须包含相同数量的观测值。如果 data 包含多个样本，则 data 内的配对观测值在样本之间独立交换。因此，在精确检验中，如果 m 是样本数，n 是每个样本中的观测数，则排列数为：

factorial(m)**n 

几种配对样本的统计检验，如威尔科克森符号秩检验和配对样本 t 检验，仅考虑两个配对元素之间的差异。因此，如果data只包含一个样本，则零假设分布是通过独立改变每个观测值的符号形成的。

警告

p 值通过计算零假设分布中与统计量观察值一样极端或更极端的元素来计算。由于使用有限精度算术，某些统计函数在理论值完全相等时返回数值上不同的值。在某些情况下，这可能导致计算 p 值时的大误差。permutation_test通过考虑与检验统计量观测值“接近”（在因子1+1e-14范围内）的零假设分布元素来防范这种情况。然而，建议用户检查零假设分布，以评估此比较方法是否合适；如果不合适，则手动计算 p 值。请参阅下面的示例。

参考文献

[1]

1. 1. Fisher. 《实验设计》，第六版（1951）。

[2]

B. Phipson 和 G. K. Smyth. “随机抽取排列 p 值不应为零：在随机绘制排列时计算精确 p 值。”《统计应用于遗传学和分子生物学》9.1（2010）。

[3]

M. D. Ernst. “排列方法：精确推断的基础”。《统计科学》（2004）。

[4]

B. Efron 和 R. J. Tibshirani. 《Bootstrap 的介绍》（1993）。

示例

假设我们希望测试两个样本是否来自同一分布。假设我们对底层分布一无所知，并且在观察数据之前，我们假设第一个样本的均值将小于第二个样本的均值。我们决定使用样本均值之差作为检验统计量，并且我们将认为 p 值为 0.05 具有统计显著性。

为了效率，我们以向量化的方式编写了定义测试统计量的函数：样本 x 和 y 可以是 ND 数组，统计量将沿着 axis 轴片段计算。

>>> import numpy as np
>>> def statistic(x, y, axis):
...     return np.mean(x, axis=axis) - np.mean(y, axis=axis) 

在收集数据后，我们计算检验统计量的观察值。

>>> from scipy.stats import norm
>>> rng = np.random.default_rng()
>>> x = norm.rvs(size=5, random_state=rng)
>>> y = norm.rvs(size=6, loc = 3, random_state=rng)
>>> statistic(x, y, 0)
-3.5411688580987266 

确实，检验统计量为负，表明 x 底层分布的真实均值小于 y 底层分布的真实均值。为了确定这种情况的概率是否由于两个样本从相同分布中抽取而偶然发生，我们执行了排列检验。

>>> from scipy.stats import permutation_test
>>> # because our statistic is vectorized, we pass `vectorized=True`
>>> # `n_resamples=np.inf` indicates that an exact test is to be performed
>>> res = permutation_test((x, y), statistic, vectorized=True,
...                        n_resamples=np.inf, alternative='less')
>>> print(res.statistic)
-3.5411688580987266
>>> print(res.pvalue)
0.004329004329004329 

在零假设下获得小于或等于观察值的检验统计量的概率为 0.4329%。这比我们选择的 5%阈值小，因此我们认为这是支持备择假设反对零假设的显著证据。

因为上述样本大小较小，permutation_test 可以执行精确检验。对于较大的样本，我们采用随机排列检验。

>>> x = norm.rvs(size=100, random_state=rng)
>>> y = norm.rvs(size=120, loc=0.3, random_state=rng)
>>> res = permutation_test((x, y), statistic, n_resamples=100000,
...                        vectorized=True, alternative='less',
...                        random_state=rng)
>>> print(res.statistic)
-0.5230459671240913
>>> print(res.pvalue)
0.00016999830001699983 

在零假设下获得小于或等于观察值的检验统计量的近似概率为 0.0225%。这同样小于我们选择的 5%阈值，因此我们再次有足够的证据来拒绝零假设，支持备择假设。

对于大样本和排列次数，结果与相应的渐近检验——独立样本 t 检验相比可比较。

>>> from scipy.stats import ttest_ind
>>> res_asymptotic = ttest_ind(x, y, alternative='less')
>>> print(res_asymptotic.pvalue)
0.00012688101537979522 

提供了进一步调查的测试统计量的排列分布。

>>> import matplotlib.pyplot as plt
>>> plt.hist(res.null_distribution, bins=50)
>>> plt.title("Permutation distribution of test statistic")
>>> plt.xlabel("Value of Statistic")
>>> plt.ylabel("Frequency")
>>> plt.show() 

如果统计量由于有限的机器精度而不准确，检查空分布至关重要。考虑以下情况：

>>> from scipy.stats import pearsonr
>>> x = [1, 2, 4, 3]
>>> y = [2, 4, 6, 8]
>>> def statistic(x, y):
...     return pearsonr(x, y).statistic
>>> res = permutation_test((x, y), statistic, vectorized=False,
...                        permutation_type='pairings',
...                        alternative='greater')
>>> r, pvalue, null = res.statistic, res.pvalue, res.null_distribution 

在这种情况下，由于数值噪声，空分布中的一些元素与检验统计量 r 的观察值不同。我们手动检查了空分布中接近检验统计量观察值的元素。

>>> r
0.8
>>> unique = np.unique(null)
>>> unique
array([-1\. , -0.8, -0.8, -0.6, -0.4, -0.2, -0.2,  0\. ,  0.2,  0.2,  0.4,
 0.6,  0.8,  0.8,  1\. ]) # may vary
>>> unique[np.isclose(r, unique)].tolist()
[0.7999999999999999, 0.8] 

如果permutation_test 在比较时过于天真，空分布中值为 0.7999999999999999 的元素将不被视为与统计量的观察值一样极端或更极端，因此计算得到的 p 值将会过小。

>>> incorrect_pvalue = np.count_nonzero(null >= r) / len(null)
>>> incorrect_pvalue
0.1111111111111111  # may vary 

相反，permutation_test 将空分布中与统计量 r 的观察值在 max(1e-14, abs(r)*1e-14) 范围内的元素视为等于 r。

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

这种比较方法预计在大多数实际情况下都是准确的，但建议用户通过检查与统计量观察值接近的空分布元素来评估此准确性。另外，考虑使用可以使用精确算术计算的统计量（例如整数统计）。

scipy.stats.bootstrap

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

scipy.stats.bootstrap(data, statistic, *, n_resamples=9999, batch=None, vectorized=None, paired=False, axis=0, confidence_level=0.95, alternative='two-sided', method='BCa', bootstrap_result=None, random_state=None)

计算统计量的双侧自举置信区间。

当method为'percentile'且alternative为'two-sided'时，根据以下过程计算自举置信区间。

1. 重新采样数据：对data中的每个样本和每个n_resamples，从原始样本中取出相同大小的随机样本（有放回）。

2. 计算统计量的自举分布：对每组重新采样计算检验统计量。

3. 确定置信区间：找到自举分布的区间，该区间为

   • 关于中位数对称且

   • 包含重新采样统计值的confidence_level。

虽然'percentile'方法最直观，但实际上很少使用。有两种更常见的方法可用，'basic'（反向百分位）和'BCa'（偏差校正和加速），它们在执行步骤 3 时有所不同。

如果data中的样本是随机从各自分布中抽取的(n)次，则bootstrap返回的置信区间将大约包含confidence_level(, \times , n)次这些分布的统计值。

参数：

data数组的序列

data的每个元素都是来自底层分布的样本。

statistic可调用函数

要计算其置信区间的统计量。statistic必须是一个可调用函数，接受len(data)个样本作为单独的参数并返回结果统计量。如果设置了vectorized为True，则statistic还必须接受一个关键字参数axis并且能够对提供的axis进行向量化计算统计量。

n_resamples整型，默认值：9999

对统计量的自举分布进行的重新采样次数。

batch整型，可选

每次对statistic进行向量化调用时处理的重新采样次数。内存使用量为 O( batch * n )，其中 n 是样本大小。默认为 None，此时 batch = n_resamples（或对于 method='BCa'，batch = max(n_resamples, n)）。

vectorized布尔型，可选

如果设置了vectorized为False，则statistic将不会传递关键字参数axis，并且预计仅计算 1D 样本的统计量。如果为True，则当传递一个 ND 样本数组时，statistic将被传递关键字参数axis，并且预计将沿着提供的axis计算统计量。如果为None（默认），则如果statistic的参数中包含axis，则vectorized将被设置为True。使用向量化统计量通常会减少计算时间。

paired布尔型，默认值：False

统计量是否将data中相应样本的元素视为配对。

axisint，默认为0

在data中样本的轴线，计算statistic的轴线。

confidence_levelfloat，默认为0.95

置信区间的置信水平。

alternative{'two-sided', ‘less’, ‘greater’}，默认为'two-sided'

选择'two-sided'（默认）用于双侧置信区间，'less'用于下限为-np.inf的单侧置信区间，'greater'用于上限为np.inf的单侧置信区间。单侧置信区间的另一边界与两侧置信区间的confidence_level两倍离 1.0 的距离相同；例如，95% 'less' 置信区间的上限与 90% 'two-sided' 置信区间的上限相同。

method{‘percentile’, ‘basic’, ‘bca’}，默认为'BCa'

是否返回'percentile'自助法置信区间（'percentile'），'basic'（也称为‘reverse’）自助法置信区间（'basic'），或修正和加速的自助法置信区间（'BCa'）。

bootstrap_resultBootstrapResult，可选

将上次调用bootstrap返回的结果对象包含在新的自助法分布中。例如，可以用来改变confidence_level，改变method，或查看执行额外重采样的效果，而不重复计算。

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

numpy.random.RandomState，可选

用于生成重采样的伪随机数生成器状态。

如果random_state为None（或np.random），则使用numpy.random.RandomState单例。如果random_state为 int，则使用种子为random_state的新的RandomState实例。如果random_state已经是Generator或RandomState实例，则使用该实例。

返回：

resBootstrapResult

一个具有属性的对象：

confidence_intervalConfidenceInterval

作为collections.namedtuple的自助法置信区间，具有low和high属性。

bootstrap_distributionndarray

自助法分布，即每个重采样的statistic值。最后一个维度对应于重采样（例如，res.bootstrap_distribution.shape[-1] == n_resamples）。

standard_errorfloat 或 ndarray

自助法标准误差，即自助法分布的样本标准偏差。

警告：

DegenerateDataWarning

当 method='BCa' 且自助法分布是退化的（例如所有元素相同）时生成。

注释

如果自助法分布是退化的（例如所有元素都相同），则置信区间的元素可能为 NaN，此时考虑使用另一 method 或检查 data，以指示其他分析可能更合适（例如所有观察结果相同）。

参考文献

[1]

B. Efron 和 R. J. Tibshirani，《自助法介绍》，Chapman & Hall/CRC，Boca Raton，FL，USA（1993）

[2]

Nathaniel E. Helwig，《自助法置信区间》，users.stat.umn.edu/~helwig/notes/bootci-Notes.pdf

[3]

自助法（统计学），维基百科，en.wikipedia.org/wiki/Bootstrapping_%28statistics%29

示例

假设我们从一个未知分布中抽取了样本数据。

>>> import numpy as np
>>> rng = np.random.default_rng()
>>> from scipy.stats import norm
>>> dist = norm(loc=2, scale=4)  # our "unknown" distribution
>>> data = dist.rvs(size=100, random_state=rng) 

我们对分布的标准偏差感兴趣。

>>> std_true = dist.std()      # the true value of the statistic
>>> print(std_true)
4.0
>>> std_sample = np.std(data)  # the sample statistic
>>> print(std_sample)
3.9460644295563863 

自助法用于近似我们期望的变异性，如果我们重复从未知分布中抽取样本并每次计算样本的统计量。它通过反复用放回地从原始样本中重新抽取值并计算每个重新抽样的统计量来实现此目的。这导致了统计量的“自助法分布”。

>>> import matplotlib.pyplot as plt
>>> from scipy.stats import bootstrap
>>> data = (data,)  # samples must be in a sequence
>>> res = bootstrap(data, np.std, confidence_level=0.9,
...                 random_state=rng)
>>> fig, ax = plt.subplots()
>>> ax.hist(res.bootstrap_distribution, bins=25)
>>> ax.set_title('Bootstrap Distribution')
>>> ax.set_xlabel('statistic value')
>>> ax.set_ylabel('frequency')
>>> plt.show() 

标准误差量化了这种变异性。它被计算为自助法分布的标准偏差。

>>> res.standard_error
0.24427002125829136
>>> res.standard_error == np.std(res.bootstrap_distribution, ddof=1)
True 

统计量的自助法分布通常近似为具有与标准误差相等的尺度的正态分布。

>>> x = np.linspace(3, 5)
>>> pdf = norm.pdf(x, loc=std_sample, scale=res.standard_error)
>>> fig, ax = plt.subplots()
>>> ax.hist(res.bootstrap_distribution, bins=25, density=True)
>>> ax.plot(x, pdf)
>>> ax.set_title('Normal Approximation of the Bootstrap Distribution')
>>> ax.set_xlabel('statistic value')
>>> ax.set_ylabel('pdf')
>>> plt.show() 

这表明，我们可以基于该正态分布的分位数构建统计量的 90%置信区间。

>>> norm.interval(0.9, loc=std_sample, scale=res.standard_error)
(3.5442759991341726, 4.3478528599786) 

由于中心极限定理，该正态近似对样本下的各种统计量和分布是准确的；然而，在所有情况下该近似并不可靠。因为 bootstrap 被设计为适用于任意的底层分布和统计量，它使用更先进的技术来生成准确的置信区间。

>>> print(res.confidence_interval)
ConfidenceInterval(low=3.57655333533867, high=4.382043696342881) 

如果我们从原始分布中抽取 1000 次样本，并为每个样本形成一个自助法置信区间，则该置信区间大约 90%的时间包含统计量的真值。

>>> n_trials = 1000
>>> ci_contains_true_std = 0
>>> for i in range(n_trials):
...    data = (dist.rvs(size=100, random_state=rng),)
...    ci = bootstrap(data, np.std, confidence_level=0.9, n_resamples=1000,
...                   random_state=rng).confidence_interval
...    if ci[0] < std_true < ci[1]:
...        ci_contains_true_std += 1
>>> print(ci_contains_true_std)
875 

我们可以一次确定所有 1000 个样本的置信区间，而不是编写循环。

>>> data = (dist.rvs(size=(n_trials, 100), random_state=rng),)
>>> res = bootstrap(data, np.std, axis=-1, confidence_level=0.9,
...                 n_resamples=1000, random_state=rng)
>>> ci_l, ci_u = res.confidence_interval 

在这里，ci_l 和 ci_u 包含 n_trials = 1000 个样本的每个置信区间。

>>> print(ci_l[995:])
[3.77729695 3.75090233 3.45829131 3.34078217 3.48072829]
>>> print(ci_u[995:])
[4.88316666 4.86924034 4.32032996 4.2822427  4.59360598] 

再次强调，约 90%的情况下包含真实值，std_true = 4。

>>> print(np.sum((ci_l < std_true) & (std_true < ci_u)))
900 

bootstrap 也可用于估计多样本统计量的置信区间，包括假设检验计算的那些。scipy.stats.mood 执行 Mood's 测试以检验等比例参数，它返回两个输出：一个统计量和一个 p 值。要获取测试统计量的置信区间，我们首先封装一个接受两个样本参数的函数，接受一个 axis 关键字参数，并仅返回统计量。

>>> from scipy.stats import mood
>>> def my_statistic(sample1, sample2, axis):
...     statistic, _ = mood(sample1, sample2, axis=-1)
...     return statistic 

这里，我们使用“百分位数”方法，默认置信水平为 95%。

>>> sample1 = norm.rvs(scale=1, size=100, random_state=rng)
>>> sample2 = norm.rvs(scale=2, size=100, random_state=rng)
>>> data = (sample1, sample2)
>>> res = bootstrap(data, my_statistic, method='basic', random_state=rng)
>>> print(mood(sample1, sample2)[0])  # element 0 is the statistic
-5.521109549096542
>>> print(res.confidence_interval)
ConfidenceInterval(low=-7.255994487314675, high=-4.016202624747605) 

标准误的 bootstrap 估计也可用。

>>> print(res.standard_error)
0.8344963846318795 

成对样本统计量也适用。例如，考虑 Pearson 相关系数。

>>> from scipy.stats import pearsonr
>>> n = 100
>>> x = np.linspace(0, 10, n)
>>> y = x + rng.uniform(size=n)
>>> print(pearsonr(x, y)[0])  # element 0 is the statistic
0.9962357936065914 

我们封装 pearsonr 函数，以便仅返回统计量。

>>> def my_statistic(x, y):
...     return pearsonr(x, y)[0] 

我们使用 paired=True 调用 bootstrap。同时，由于 my_statistic 未矢量化以计算给定轴上的统计量，我们传入 vectorized=False。

>>> res = bootstrap((x, y), my_statistic, vectorized=False, paired=True,
...                 random_state=rng)
>>> print(res.confidence_interval)
ConfidenceInterval(low=0.9950085825848624, high=0.9971212407917498) 

结果对象可以传回 bootstrap 进行额外的重采样：

>>> len(res.bootstrap_distribution)
9999
>>> res = bootstrap((x, y), my_statistic, vectorized=False, paired=True,
...                 n_resamples=1001, random_state=rng,
...                 bootstrap_result=res)
>>> len(res.bootstrap_distribution)
11000 

或更改置信区间选项：

>>> res2 = bootstrap((x, y), my_statistic, vectorized=False, paired=True,
...                  n_resamples=0, random_state=rng, bootstrap_result=res,
...                  method='percentile', confidence_level=0.9)
>>> np.testing.assert_equal(res2.bootstrap_distribution,
...                         res.bootstrap_distribution)
>>> res.confidence_interval
ConfidenceInterval(low=0.9950035351407804, high=0.9971170323404578) 

无需重复计算原始 bootstrap 分布。

scipy.stats.MonteCarloMethod

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

class scipy.stats.MonteCarloMethod(n_resamples=9999, batch=None, rvs=None)

用于蒙特卡洛假设检验的配置信息。

可将此类的实例传递给某些假设检验函数的method参数，以执行假设检验的蒙特卡洛版本。

属性：

n_resamples整数，可选

要抽取的蒙特卡洛样本数。默认值为 9999。

batch整数，可选

在每次对统计量进行向量化调用时要处理的蒙特卡洛样本数。当统计量被向量化时，批量大小 >>1 通常更快，但内存使用量与批量大小呈线性关系。默认值为None，将所有样本在单个批次中处理。

rvs可调用对象或者可调用对象的元组，可选

一个可调用或者一系列在零假设下生成随机变量的可调用对象。每个rvs的元素必须是一个接受关键字参数size（例如rvs(size=(m, n))）并返回该形状的 N 维数组样本的可调用对象。如果rvs是一个序列，则rvs中的可调用对象数量必须与在使用MonteCarloMethod的假设检验中传递给样本数相匹配。默认值为None，此时假设检验函数选择值以匹配假设检验的标准版本。例如，scipy.stats.pearsonr的零假设通常是样本是从标准正态分布中抽取的，因此rvs = (rng.normal, rng.normal)，其中rng = np.random.default_rng()。

scipy.stats.PermutationMethod

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

class scipy.stats.PermutationMethod(n_resamples=9999, batch=None, random_state=None)

一种置换假设检验的配置信息。

此类的实例可以传递到某些假设检验函数的method参数中，以执行假设检验的置换版本。

属性：

n_resamplesint，可选

要执行的重采样次数。默认值为 9999。

batchint，可选

在每次向量化调用统计量时处理的重采样次数。当统计量被向量化时，批处理大小 >>1 通常更快，但内存使用量与批处理大小线性扩展。默认为None，即在单个批处理中处理所有重采样。

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

numpy.random.RandomState，可选

用于生成重采样的伪随机数生成器状态。

如果random_state已经是Generator或RandomState实例，则使用该实例。如果random_state是一个整数，则使用一个新的RandomState实例，并使用random_state进行种子化。如果random_state为None（默认），则使用numpy.random.RandomState单例。

scipy.stats.BootstrapMethod

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

class scipy.stats.BootstrapMethod(n_resamples=9999, batch=None, random_state=None, method='BCa')

自举置信区间的配置信息。

此类的实例可以传递到某些置信区间方法的method参数中，以生成自举置信区间。

属性：

n_resamplesint, optional

要执行的重采样次数。默认为 9999。

batchint, optional

在每个矢量化调用统计量中处理的重采样次数。当统计量被矢量化时，批量大小>>1 通常更快，但内存使用量与批量大小成线性关系。默认为None，即在单个批次中处理所有重采样。

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

numpy.random.RandomState}, optional

用于生成重采样的伪随机数生成器状态。

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

method{‘bca’, ‘percentile’, ‘basic’}

是否使用‘percentile’自举法（‘percentile’），‘basic’（又名‘reverse’）自举法（‘basic’），或者校正和加速的自举法（‘BCa’，默认值）。

scipy.stats.combine_pvalues

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

scipy.stats.combine_pvalues(pvalues, method='fisher', weights=None)

结合对同一假设相关的独立测试的 p 值。

这些方法仅用于组合基于连续分布的假设检验的 p 值。

每种方法假设在零假设下，p 值是独立采样且均匀分布于区间 [0, 1]。计算一个检验统计量（每种方法不同），并根据此检验统计量在零假设下的分布计算组合 p 值。

参数：

pvaluesarray_like，1-D

假设来自基于连续分布的独立测试的 p 值数组。

方法{‘fisher’, ‘pearson’, ‘tippett’, ‘stouffer’, ‘mudholkar_george’}

要使用的方法名称来组合 p 值。

可用的方法有（详见笔记）：

• ‘fisher’：费舍尔方法（费舍尔组合概率检验）

• ‘pearson’：皮尔逊方法

• ‘mudholkar_george’：穆德霍尔卡和乔治的方法

• ‘tippett’：提普特方法

• ‘stouffer’：斯托弗的 Z 分数方法

weightsarray_like，1-D，可选

仅用于斯托弗的 Z 分数方法的权重数组。

返回：

res显著性结果

一个包含属性的对象：

statisticfloat

指定方法计算的统计量。

pvaluefloat

组合 p 值。

笔记

如果此函数应用于具有离散统计量（例如任何秩测试或列联表测试）的测试，它将产生系统性错误的结果，例如费舍尔方法将系统性高估 p 值[1]。对于大样本量时，离散分布近似连续时，这个问题变得不那么严重。

方法之间的差异可以通过它们的统计量和在考虑显著性时强调 p 值的哪些方面来最好地说明[2]。例如，强调大的 p 值的方法对强假阴性和真阴性更为敏感；相反，侧重于小的 p 值的方法对阳性敏感。

• 费舍尔方法的统计量（也称为费舍尔组合概率检验）[3] 是 (-2\sum_i \log(p_i))，它等价于（作为一个检验统计量）各个 p 值的乘积：(\prod_i p_i)。在零假设下，这一统计量服从 (\chi²) 分布。此方法强调小的 p 值。

• 皮尔逊方法使用 (-2\sum_i\log(1-p_i))，它等价于 (\prod_i \frac{1}{1-p_i}) [2]。因此，它强调大的 p 值。

• Mudholkar 和 George 通过平均他们的统计方法在 Fisher 和 Pearson 方法之间做出妥协[4]。他们的方法强调极端的 p 值，无论是接近 1 还是 0。

• Stouffer 方法[5]使用 Z 分数和统计量：(\sum_i \Phi^{-1} (p_i))，其中(\Phi)是标准正态分布的累积分布函数。该方法的优势在于可以简单地引入权重，这可以使 Stouffer 方法在来自不同大小研究的 p 值时比 Fisher 方法更有效[6] [7]。

• Tippett 方法使用最小的 p 值作为统计量。（请注意，这个最小值不是组合 p 值。）

Fisher 方法可以扩展到组合来自相关测试的 p 值[8]。目前未实现的扩展方法包括 Brown 方法和 Kost 方法。

新版本 0.15.0 中的新内容。

参考文献

[1]

Kincaid, W. M., “The Combination of Tests Based on Discrete Distributions.” Journal of the American Statistical Association 57, no. 297 (1962), 10-19.

[2] (1,2)

Heard, N. and Rubin-Delanchey, P. “Choosing between methods of combining p-values.” Biometrika 105.1 (2018): 239-246.

[3]

en.wikipedia.org/wiki/Fisher%27s_method

[4]

George, E. O., and G. S. Mudholkar. “On the convolution of logistic random variables.” Metrika 30.1 (1983): 1-13.

[5]

en.wikipedia.org/wiki/Fisher%27s_method#Relation_to_Stouffer.27s_Z-score_method

[6]

Whitlock, M. C. “Combining probability from independent tests: the weighted Z-method is superior to Fisher’s approach.” Journal of Evolutionary Biology 18, no. 5 (2005): 1368-1373.

[7]

Zaykin, Dmitri V. “Optimally weighted Z-test is a powerful method for combining probabilities in meta-analysis.” Journal of Evolutionary Biology 24, no. 8 (2011): 1836-1841.

[8]

en.wikipedia.org/wiki/Extensions_of_Fisher%27s_method

例子

假设我们希望使用 Fisher 方法（默认）来组合相同零假设的四个独立测试的 p 值。

>>> from scipy.stats import combine_pvalues
>>> pvalues = [0.1, 0.05, 0.02, 0.3]
>>> combine_pvalues(pvalues)
SignificanceResult(statistic=20.828626352604235, pvalue=0.007616871850449092) 

当各个 p 值具有不同的权重时，请考虑 Stouffer 方法。

>>> weights = [1, 2, 3, 4]
>>> res = combine_pvalues(pvalues, method='stouffer', weights=weights)
>>> res.pvalue
0.009578891494533616 

scipy.stats.false_discovery_control

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

scipy.stats.false_discovery_control(ps, *, axis=0, method='bh')

调整 p 值以控制假发现率。

假发现率（FDR）是被拒绝的空假设中实际为真的比例的期望值。如果在调整后 p 值低于指定水平时拒绝空假设，则假发现率在该水平上得到控制。

参数：

ps：1D array_like

需要调整的 p 值。元素必须是介于 0 和 1 之间的实数。

axis：int

执行调整的轴。沿每个轴切片独立执行调整。如果 axis 为 None，则在执行调整之前对 ps 进行展平。

method：{‘bh’，‘by’}

应用的假发现率控制程序：'bh'指的是本雅明-霍克伯格[1]（方程 1），'by'指的是本雅明-耶库提耶里[2]（定理 1.3）。后者更为保守，但确保即使 p 值不是来自独立测试，也能控制假发现率。

返回：

ps_adjusted：array_like

调整后的 p 值。如果这些值低于指定水平时拒绝空假设，则假发现率在该水平上得到控制。

参见

combine_pvalues

statsmodels.stats.multitest.multipletests

注：

在多重假设检验中，假发现控制程序往往比家族误差率控制程序（例如 Bonferroni 校正[1]）提供更高的功效。

如果 p 值对应于独立测试（或具有“正回归依赖性”的测试[2]），拒绝 Benjamini-Hochberg 调整后 p 值低于 (q) 的空假设可以控制假发现率在小于或等于 (q m_0 / m) 的水平上，其中 (m_0) 是真空假设的数量，(m) 是测试的总空假设数量。即使对于依赖测试，当根据更保守的 Benjaminini-Yekutieli 程序进行调整时，情况也是如此。

本函数生成的调整后的 p 值可与 R 函数 p.adjust 和 statsmodels 函数 statsmodels.stats.multitest.multipletests 生成的相比较。请考虑后者以获取更高级的多重比较校正方法。

参考文献

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

Benjamini, Yoav, 和 Yosef Hochberg. “控制假发现率：多重检验的实用和强大方法。” 王立统计学会系列 B (方法论) 57.1 (1995): 289-300.

[2] (1,2)

Benjamini, Yoav, 和 Daniel Yekutieli. “控制相关性下的多重检验假阳率。” 统计学年鉴 (2001): 1165-1188.

[3]

TileStats. FDR - Benjamini-Hochberg explained - Youtube. www.youtube.com/watch?v=rZKa4tW2NKs.

[4]

Neuhaus, Karl-Ludwig, 等。“rt-PA-APSAC 通透性研究(TAPS)：急性心肌梗死中通过 rt-PA 前负荷治疗改善溶栓治疗效果。” 美国心脏病学会杂志 19.5 (1992): 885-891.

例子

我们遵循[1]的例子。

在心肌梗死中，利用重组组织型纤溶酶原激活剂(rt-PA)和苯乙酰化的纤溶酶原激活剂(APSAC)的溶栓治疗已被证明能够降低死亡率。[4]在一项随机多中心试验中，研究了新的 rt-PA 前负荷治疗与标准 APSAC 方案治疗在 421 例急性心肌梗死患者中的效果。

研究中测试了四个假设家族，最后一个是“心脏和其他事件在溶栓治疗开始后”。在这个假设家族中，可能需要 FDR 控制，因为如果前负荷治疗仅与先前治疗相当，则不宜得出前者更佳的结论。

此家族中 15 个假设对应的 p 值如下：

>>> ps = [0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344,
...       0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000] 

如果所选显著性水平为 0.05，我们可能会倾向于拒绝前九个 p 值对应的零假设，因为前九个 p 值低于所选显著性水平。然而，这会忽略“多重性”的问题：如果我们未能纠正多重比较的事实，我们更有可能错误地拒绝真实的零假设。

解决多重性问题的一种方法是控制家族错误率(FWER)，即在零假设实际为真时拒绝的比率。这种类型的常见程序是 Bonferroni 校正[1]。我们首先将 p 值乘以测试的假设数。

>>> import numpy as np
>>> np.array(ps) * len(ps)
array([1.5000e-03, 6.0000e-03, 2.8500e-02, 1.4250e-01, 3.0150e-01,
 4.1700e-01, 4.4700e-01, 5.1600e-01, 6.8850e-01, 4.8600e+00,
 6.3930e+00, 8.5785e+00, 9.7920e+00, 1.1385e+01, 1.5000e+01]) 

控制 FWER 在 5%水平下，我们仅拒绝调整后的 p 值小于 0.05 的假设。在这种情况下，只有与前三个 p 值相关的假设可以被拒绝。根据[1]，这三个假设涉及“过敏反应”和“出血的两个不同方面”。

另一种方法是控制虚假发现率：预期被拒绝的零假设中实际为真的比例。这种方法的优势在于，通常提供更大的功效：在零假设确实为假时，拒绝零假设的增加率。为了将虚假发现率控制在 5%以内，我们采用 Benjamini-Hochberg p 值调整方法。

>>> from scipy import stats
>>> stats.false_discovery_control(ps)
array([0.0015    , 0.003     , 0.0095    , 0.035625  , 0.0603    ,
 0.06385714, 0.06385714, 0.0645    , 0.0765    , 0.486     ,
 0.58118182, 0.714875  , 0.75323077, 0.81321429, 1\.        ]) 

现在，四 个调整后的 p 值第一次低于 0.05，因此我们将拒绝与这些 四 个 p 值对应的零假设。特别重要的是第四个零假设的拒绝，因为它导致了结论：新治疗方法的“住院死亡率显著降低”。

scipy.stats.boxcox

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

scipy.stats.boxcox(x, lmbda=None, alpha=None, optimizer=None)

返回通过 Box-Cox 幂变换转换的数据集。

参数：

xndarray

要转换的输入数组。

如果lmbda不是 None，则这是scipy.special.boxcox的别名。如果x < 0，返回 nan；如果x == 0且lmbda < 0，返回-inf。

如果lmbda为 None，则数组必须是正的、一维的且非常数。

lmbdascalar，可选

如果lmbda为 None（默认），则找到最大化对数似然函数的lmbda值并将其作为第二个输出参数返回。

如果lmbda不是 None，则对该值进行转换。

alphafloat，可选

如果lmbda为 None 且alpha不为 None（默认），则将lmbda的100 * (1-alpha)%置信区间作为第三个输出参数返回。必须介于 0.0 和 1.0 之间。

如果lmbda不是 None，将忽略alpha。

optimizercallable，可选

如果lmbda为 None，则optimizer是用于找到最小化负对数似然函数的lmbda值的标量优化器。optimizer是一个接受一个参数的可调用对象：

funcallable

目标函数，用于在提供的lmbda值处评估负对数似然函数。

并返回一个对象，例如scipy.optimize.OptimizeResult，其中在属性x中保存了最优的lmbda值。

更多信息请参见boxcox_normmax中的示例或scipy.optimize.minimize_scalar的文档。

如果lmbda不是 None，则忽略optimizer。

返回：

boxcoxndarray

Box-Cox 幂变换的数组。

maxlogfloat，可选

如果lmbda参数为 None，则第二个返回参数是最大化对数似然函数的lmbda值。

**(min_ci, max_ci)**float 元组，可选

如果lmbda参数为 None 且alpha不为 None，则返回的这个浮点数元组表示给定alpha的最小和最大置信限。

另请参见

probplot，boxcox_normplot，boxcox_normmax，boxcox_llf

注意

Box-Cox 变换由以下提供：

y = (x**lmbda - 1) / lmbda,  for lmbda != 0
    log(x),                  for lmbda = 0 

boxcox 要求输入数据为正数。有时 Box-Cox 变换提供一个移动参数以实现此目的；boxcox 并不提供此类移动参数。这样的移动参数等同于在调用 boxcox 之前向 x 添加一个正常数。

当提供 alpha 时返回的置信限给出了以下区间：

[llf(\hat{\lambda}) - llf(\lambda) < \frac{1}{2}\chi²(1 - \alpha, 1),]

这里的 llf 表示对数似然函数，(\chi²) 表示卡方函数。

参考文献

G.E.P. Box 和 D.R. Cox，《转换的分析》，《皇家统计学会 B》杂志，26，211-252（1964）。

示例

>>> from scipy import stats
>>> import matplotlib.pyplot as plt 

我们从非正态分布生成一些随机变量，并制作概率图来展示其在尾部的非正态性：

>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(211)
>>> x = stats.loggamma.rvs(5, size=500) + 5
>>> prob = stats.probplot(x, dist=stats.norm, plot=ax1)
>>> ax1.set_xlabel('')
>>> ax1.set_title('Probplot against normal distribution') 

现在我们使用 boxcox 对数据进行转换，使其尽可能接近正态分布：

>>> ax2 = fig.add_subplot(212)
>>> xt, _ = stats.boxcox(x)
>>> prob = stats.probplot(xt, dist=stats.norm, plot=ax2)
>>> ax2.set_title('Probplot after Box-Cox transformation') 

>>> plt.show()