SciPy-1-12-中文文档-二十七-
SciPy 1.12 中文文档（二十七）
scipy.stats.dunnett
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.dunnett.html#scipy.stats.dunnett
scipy.stats.dunnett(*samples, control, alternative='two-sided', random_state=None)
Dunnett’s 测试：多个组的平均值与控制组的比较。
这是 Dunnett 原始的单步测试实现，如[1]所述。
参数：
sample1, sample2, …1D 数组样本
每个实验组的样本测量。
control1D 数组样本
控制组的样本测量。
alternative，可选
定义备择假设。
零假设是样本分布和控制组分布的平均值相等。可用以下备择假设（默认为‘双边’）：
-
‘two-sided’: 样本和控制组的分布平均值不等。
-
‘less’: 样本分布的平均值小于控制组分布的平均值。
-
‘greater’: 样本分布的平均值大于控制组分布的平均值。
random_state{None, int, numpy.random.Generator}, 可选
如果random_state是 int 或 None，则使用np.random.default_rng(random_state)创建一个新的numpy.random.Generator。如果random_state已经是Generator实例，则使用提供的实例。
随机数生成器用于控制多元 t 分布的随机化拟蒙特卡罗积分。
返回：
包含属性的对象：
统计量的浮点数 ndarray
每次比较的测试计算统计量。索引i处的元素是组i与控制组之间的统计量。
p 值的浮点数 ndarray
每次比较的测试的计算 p 值。索引i处的元素是组i与控制组之间的 p 值。
以及以下方法：
confidence_interval(confidence_level=0.95)：
计算组的平均值与控制组加减允许范围的差异。
另见
执行平均值的成对比较。
注意事项
像独立样本 t 检验一样，邓特氏检验[1]用于对样本抽样分布的均值进行推断。然而，当以固定显著性水平执行多个 t 检验时，“家族误差率” - 在至少一个测试中错误拒绝零假设的概率 - 将超过显著性水平。邓特氏检验旨在在控制家族误差率的同时进行多重比较。
邓特氏检验比较多个实验组与单一对照组的均值。Tukey 的 Honestly Significant Difference Test 是另一种控制家族误差率的多重比较测试，但tukey_hsd进行所有成对组间比较。当不需要实验组间的成对比较时，邓特氏检验由于具有更高的功效而更可取。
此测试的使用依赖于几个假设。
-
观察在组内和组间是独立的。
-
每组观察值符合正态分布。
-
抽样分布具有相同的有限方差。
参考资料
[1] (1,2,3)
Charles W. Dunnett. “A Multiple Comparison Procedure for Comparing Several Treatments with a Control.” 美国统计协会杂志, 50:272, 1096-1121, DOI:10.1080/01621459.1955.10501294, 1955.
例子
在[1]中，研究了药物对三组动物血细胞计数测量的影响。
下表总结了实验结果，其中两组接受不同药物，而一组作为对照。记录了血细胞计数（每立方毫米百万细胞数）：
>>> import numpy as np
>>> control = np.array([7.40, 8.50, 7.20, 8.24, 9.84, 8.32])
>>> drug_a = np.array([9.76, 8.80, 7.68, 9.36])
>>> drug_b = np.array([12.80, 9.68, 12.16, 9.20, 10.55])
我们想看看各组的平均值是否显著不同。首先，通过箱线图进行视觉检查。
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(1, 1)
>>> ax.boxplot([control, drug_a, drug_b])
>>> ax.set_xticklabels(["Control", "Drug A", "Drug B"])
>>> ax.set_ylabel("mean")
>>> plt.show()
注意药物 A 组和对照组的重叠四分位范围以及药物 B 组和对照组之间的明显分离。
接下来，我们将使用邓特氏检验来评估组间均值差异是否显著，同时控制家族误差率：即可能发生任何虚假发现的概率。设定零假设为实验组与对照组均值相同，备择假设为实验组与对照组均值不同。我们将考虑 5%的家族误差率是可接受的，因此我们选择 0.05 作为显著性阈值。
>>> from scipy.stats import dunnett
>>> res = dunnett(drug_a, drug_b, control=control)
>>> res.pvalue
array([0.62004941, 0.0059035 ]) # may vary
在组 A 和对照组之间进行比较的 p 值超过了 0.05，因此我们在这个比较中不拒绝原假设。然而，在组 B 和对照组之间进行比较的 p 值小于 0.05，因此我们认为实验结果支持备择假设：组 B 的均值与对照组不同。
scipy.stats.kruskal
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.kruskal.html#scipy.stats.kruskal
scipy.stats.kruskal(*samples, nan_policy='propagate', axis=0, keepdims=False)
计算独立样本的 Kruskal-Wallis H 检验。
Kruskal-Wallis H 检验检验假设，即所有组的总体中位数相等。这是方差分析的非参数版本。该检验适用于两个或更多个独立样本，这些样本可能具有不同的大小。请注意，拒绝原假设并不表示哪些组之间不同。需要进行事后比较以确定哪些组不同。
参数：
sample1, sample2, …array_like
可以将两个或更多个数组与样本测量值作为参数给出。样本必须是一维的。
nan_policy
定义如何处理输入的 NaN。
-
propagate: 如果在计算统计量的轴切片（例如行）中存在 NaN，则输出的相应条目将为 NaN。 -
omit: 在执行计算时将忽略 NaN。如果沿计算统计量的轴切片中剩余的数据不足，则输出的相应条目将为 NaN。 -
raise: 如果存在 NaN，则会引发ValueError。
axisint 或 None，默认值：0
如果是整数，则是计算统计量的输入轴。输入的每个轴切片（例如行）的统计量将出现在输出的相应元素中。如果为None，则在计算统计量之前将对输入进行展平。
keepdimsbool，默认值：False
如果设置为 True，则减少的轴将作为尺寸为一的维度留在结果中。使用此选项，结果将正确传播到输入数组。
返回：
statisticfloat
经校正的 Kruskal-Wallis H 统计量，考虑到并列。
pvaluefloat
使用假设 H 服从卡方分布的测试的 p 值。返回的 p 值是在 H 处评估的卡方分布的生存函数。
另请参见
f_oneway
1-way ANOVA.
mannwhitneyu
两个样本的曼-惠特尼秩和检验。
friedmanchisquare
重复测量的弗里德曼检验。
注意事项
由于假设 H 服从卡方分布，每个组中的样本数量不能太少。典型规则是每个样本必须至少有 5 次测量。
从 SciPy 1.9 开始，np.matrix 输入（不建议在新代码中使用）在执行计算之前会转换为 np.ndarray。在这种情况下，输出将是适当形状的标量或 np.ndarray，而不是 2D 的 np.matrix。类似地，尽管掩码数组的掩码元素被忽略，输出将是标量或 np.ndarray，而不是具有 mask=False 的掩码数组。
参考
[1]
W. H. Kruskal & W. W. Wallis，《使用排名进行单因素方差分析》，《美国统计协会杂志》，第 47 卷，260 期，第 583-621 页，1952 年。
[2]
en.wikipedia.org/wiki/Kruskal-Wallis_one-way_analysis_of_variance
示例
>>> from scipy import stats
>>> x = [1, 3, 5, 7, 9]
>>> y = [2, 4, 6, 8, 10]
>>> stats.kruskal(x, y)
KruskalResult(statistic=0.2727272727272734, pvalue=0.6015081344405895)
>>> x = [1, 1, 1]
>>> y = [2, 2, 2]
>>> z = [2, 2]
>>> stats.kruskal(x, y, z)
KruskalResult(statistic=7.0, pvalue=0.0301973834223185)
scipy.stats.alexandergovern
scipy.stats.alexandergovern(*samples, nan_policy='propagate')
执行 Alexander Govern 检验。
Alexander-Govern 近似检验在方差异质性情况下测试 k 个独立均值的相等性。该检验适用于来自两个或多个组的样本，可能具有不同的大小。
参数：
sample1, sample2, …array_like
每组的样本测量。至少需要两个样本。
nan_policy，可选
定义如何处理包含 NaN 的输入。提供以下选项（默认为‘propagate’）：
-
‘propagate’: 返回 NaN
-
‘raise’: 抛出错误
-
‘omit’: 在计算时忽略 NaN 值
返回：
resAlexanderGovernResult
具有属性的对象：
statisticfloat
测试的计算 A 统计量。
pvaluefloat
从卡方分布中得到的关联 p 值。
警告：
ConstantInputWarning
如果输入是常数数组，则引发错误。在这种情况下，统计量未定义，因此返回np.nan。
另请参阅
f_oneway
单因素方差分析
注意事项
此检验的使用依赖于几个假设。
-
样本是独立的。
-
每个样本来自正态分布的总体。
-
与
f_oneway不同，此检验不假设同方差性，而是放宽了等方差性的假设。
输入样本必须是有限的、一维的，并且大小大于一。
参考文献
[1]
Alexander, Ralph A., 和 Diane M. Govern. “A New and Simpler Approximation for ANOVA under Variance Heterogeneity.” Journal of Educational Statistics, vol. 19, no. 2, 1994, pp. 91-101. JSTOR, www.jstor.org/stable/1165140. 访问日期：2020 年 9 月 12 日。
示例
>>> from scipy.stats import alexandergovern
这里提供了来自美国四个城市九家最大银行新车贷款年利率的一些数据，取自国家标准技术研究所的 ANOVA 数据集。
我们使用alexandergovern来检验所有城市的平均年利率百分比是否相同的零假设，对立假设是不是所有城市的平均年利率百分比都不相同。我们决定在显著性水平为 5%时拒绝零假设，支持备择假设。
>>> atlanta = [13.75, 13.75, 13.5, 13.5, 13.0, 13.0, 13.0, 12.75, 12.5]
>>> chicago = [14.25, 13.0, 12.75, 12.5, 12.5, 12.4, 12.3, 11.9, 11.9]
>>> houston = [14.0, 14.0, 13.51, 13.5, 13.5, 13.25, 13.0, 12.5, 12.5]
>>> memphis = [15.0, 14.0, 13.75, 13.59, 13.25, 12.97, 12.5, 12.25,
... 11.89]
>>> alexandergovern(atlanta, chicago, houston, memphis)
AlexanderGovernResult(statistic=4.65087071883494,
pvalue=0.19922132490385214)
p 值为 0.1992，表示在零假设下观察到这样一个极端值的几率接近 20%。这超过了 5%，因此我们不拒绝零假设，支持备择假设。
scipy.stats.fligner
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.fligner.html#scipy.stats.fligner
scipy.stats.fligner(*samples, center='median', proportiontocut=0.05, axis=0, nan_policy='propagate', keepdims=False)
执行 Fligner-Killeen 检验以检验方差的相等性。
Fligner 的检验检验的是所有输入样本来自方差相等的总体的零假设。当总体相同时，Fligner-Killeen 的检验是分布自由的 [2]。
参数：
sample1, sample2, …array_like
样本数据数组。不需要具有相同的长度。
center，可选
控制在计算检验统计量时使用数据的函数的关键字参数。默认为 ‘median’。
proportiontocutfloat，可选
当 center 为 ‘trimmed’ 时，这指定要从每端剪切的数据点的比例（参见 scipy.stats.trim_mean）。默认值为 0.05。
axisint 或 None，默认值：0
如果是整数，则沿着计算统计量的输入轴（例如行）的轴。输入的每个轴切片的统计量将出现在输出的相应元素中。如果为 None，则在计算统计量之前将展平输入。
nan_policy
定义如何处理输入的 NaN 值。
-
propagate：如果在计算统计量的轴切片（例如行）中存在 NaN，则输出的相应条目将是 NaN。 -
omit：在执行计算时将省略 NaN。如果在计算统计量的轴切片中剩余的数据不足，则输出的相应条目将是 NaN。 -
raise：如果存在 NaN，则会引发ValueError。
keepdimsbool，默认值：False
如果设置为 True，则减少的轴将作为大小为一的维度保留在结果中。使用此选项，结果将正确地与输入数组广播。
返回：
statisticfloat
检验统计量。
pvaluefloat
假设检验的 p 值。
参见
正态样本中 k 方差相等性的参数检验
一个关于 k 方差相等性的鲁棒参数检验
注意
与 Levene 的检验类似，Fligner 的检验有三个变体，它们在测试中使用的集中趋势测量方法不同。详见 levene 获取更多信息。
康诺弗等人（1981）通过广泛的模拟研究了许多现有的参数和非参数检验方法，他们得出结论，Fligner 和 Killeen（1976）以及 Levene（1960）提出的检验方法在对正态性偏差和功效的鲁棒性方面似乎更为优越 [3]。
从 SciPy 1.9 开始，不推荐使用的 np.matrix 输入在进行计算之前将被转换为 np.ndarray。在这种情况下，输出将是适当形状的标量或 np.ndarray，而不是二维 np.matrix。同样，忽略掩码数组的掩码元素后，输出将是一个标量或 np.ndarray，而不是具有 mask=False 的掩码数组。
参考文献
[1]
Park, C. 和 Lindsay, B. G. (1999). 鲁棒的尺度估计和基于二次推理函数的假设检验。宾夕法尼亚州立大学可能性研究中心技术报告 #99-03。cecas.clemson.edu/~cspark/cv/paper/qif/draftqif2.pdf
[2]
Fligner, M.A. 和 Killeen, T.J. (1976). Distribution-free two-sample tests for scale. ‘Journal of the American Statistical Association.’ 71(353), 210-213.
[3]
Park, C. 和 Lindsay, B. G. (1999). 鲁棒的尺度估计和基于二次推理函数的假设检验。宾夕法尼亚州立大学可能性研究中心技术报告 #99-03。
[4]
Conover, W. J., Johnson, M. E. 和 Johnson M. M. (1981). 各种方差同质性检验的比较研究，适用于外大陆架招标数据。技术统计学，23(4)，351-361。
[5]
C.I. BLISS (1952)，生物测定的统计学：特别参考维生素，第 499-503 页，DOI:10.1016/C2013-0-12584-6。
[6]
B. Phipson 和 G. K. Smyth. “当置换是随机抽取时，永远不应为零：计算精确的 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
... ])
fligner 统计量对样本之间方差的差异很敏感。
>>> from scipy import stats
>>> res = stats.fligner(small_dose, medium_dose, large_dose)
>>> res.statistic
1.3878943408857916
当方差存在显著差异时，统计量的值往往较高。
我们可以通过比较统计量的观察值与零假设的分布来测试组间方差的不等性：在零假设下，三组的总体方差相等的统计值的分布。
对于这个检验，零假设分布如下，遵循卡方分布。
>>> import matplotlib.pyplot as plt
>>> k = 3 # number of samples
>>> dist = stats.chi2(df=k-1)
>>> val = np.linspace(0, 8, 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("Fligner Test Null Distribution")
... ax.set_xlabel("statistic")
... ax.set_ylabel("probability density")
... ax.set_xlim(0, 8)
... ax.set_ylim(0, 0.5)
>>> plot(ax)
>>> plt.show()
比较由 p 值量化：在零分布中大于或等于统计量的观察值的比例。
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> pvalue = 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, (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.49960016501182125
如果 p 值“小” - 也就是说，从具有相同方差的分布中抽取数据并产生统计量的极端值的概率较低 - 这可能被视为反对零假设的证据，支持备择假设：这些组的方差不相等。注意：
-
反之不成立；也就是说，该检验不能用来支持零假设。
-
被视为“小”的值的阈值是在数据分析之前应该做出的选择，考虑到假阳性（错误拒绝零假设）和假阴性（未能拒绝错误的零假设）的风险[6]。
-
小的 p 值并不是支持 大 效应的证据；相反，它们只能提供“显著”效应的证据，意味着在零假设下发生这些情况的可能性很小。
注意，卡方分布提供了零分布的渐近近似。对于小样本，执行置换检验可能更合适：在零假设下，三个样本都是从同一总体抽取的，每个测量在三个样本中的观察概率相同。因此，我们可以通过在观察数据的许多随机分区下计算统计量来形成随机化的零分布。
>>> def statistic(*samples):
... return stats.fligner(*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, 8, 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.4332 # may vary
请注意，这里计算的 p 值与上述 fligner 返回的渐近近似存在显著分歧。从置换检验中可以严谨地推断的统计推断是有限的；尽管如此，在许多情况下，它们可能是首选的方法[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.fligner(a, b, c)
>>> p
0.00450826080004775
小的 p 值表明这些群体的方差不相等。
鉴于 b 的样本方差远大于 a 和 c 的样本方差，这并不令人意外：
>>> [np.var(x, ddof=1) for x in [a, b, c]]
[0.007054444444444413, 0.13073888888888888, 0.008890000000000002]
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, optional
在测试中使用数据的哪个函数。默认为 'median'。
proportiontocutfloat, optional
当 center 为 ‘trimmed’ 时，这给出了要从每端裁剪的数据点的比例。（见 scipy.stats.trim_mean.）默认为 0.05。
axisint or None, default: 0
如果是整数，则是沿其计算统计量的输入的轴。输入的每个轴切片（例如行）的统计量将出现在输出的相应元素中。如果 None，则在计算统计量之前将展平输入。
nan_policy
定义如何处理输入的 NaN。
-
propagate: 如果轴（例如行）上存在 NaN，则计算统计量时对应的输出条目将为 NaN。 -
omit: 在执行计算时将省略 NaN。如果在计算统计量时轴切片上剩余的数据不足，则对应的输出条目将为 NaN。 -
raise: 如果存在 NaN，则会引发ValueError。
keepdimsbool, default: False
如果设置为 True，则被减少的轴将作为具有大小为一的维度保留在结果中。使用此选项，结果将正确地广播到输入数组。
返回：
statisticfloat
检验统计量。
pvaluefloat
测试的 p 值。
另请参阅
正态样本中 k 个方差的非参数检验
正态样本中 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
定义如何处理输入 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
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，可选
定义如何处理输入包含 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
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
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
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
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，可选
要执行的置换类型，符合零假设的要求。前两种置换类型适用于配对样本统计量，其中所有样本包含相同数量的观测值，并且沿着 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，可选
用于计算 p 值的备择假设。对于每个备择假设，p 值的定义如下。
-
'greater'：空值分布中大于或等于测试统计量观察值的百分比。 -
'less'：空值分布中小于或等于测试统计量观察值的百分比。 -
'two-sided'（默认）：上述 p 值之一的两倍较小的值。
注意，随机化测试的 p 值是根据[2]和[3]中建议的保守（过估计）近似计算的，而不是建议的无偏估计器[4]。也就是说，在计算随机化空值分布中与测试统计量观察值一样极端的比例时，分子和分母的值都增加了一。这种调整的解释是，测试统计量的观察值总是作为随机化空值分布的一个元素。用于双边 p 值的约定不是普遍适用的；如果喜欢不同的定义，则返回观察到的测试统计量和空值分布。
axisint，默认值：0
(广播)样本的轴，用于计算统计量。如果样本具有不同维数，则在考虑 axis 之前，对具有较少维度的样本前置单例维度。
random_state{None, int, numpy.random.Generator，
用于生成排列的伪随机数生成器状态。
如果 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]
-
- 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'时，根据以下过程计算自举置信区间。
-
重新采样数据：对data中的每个样本和每个n_resamples，从原始样本中取出相同大小的随机样本（有放回）。
-
计算统计量的自举分布：对每组重新采样计算检验统计量。
-
确定置信区间：找到自举分布的区间，该区间为
-
关于中位数对称且
-
包含重新采样统计值的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，默认为'BCa'
是否返回'percentile'自助法置信区间（'percentile'），'basic'（也称为‘reverse’）自助法置信区间（'basic'），或修正和加速的自助法置信区间（'BCa'）。
bootstrap_resultBootstrapResult，可选
将上次调用bootstrap返回的结果对象包含在新的自助法分布中。例如，可以用来改变confidence_level，改变method，或查看执行额外重采样的效果，而不重复计算。
random_state{None, int, numpy.random.Generator,
用于生成重采样的伪随机数生成器状态。
如果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
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
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,
用于生成重采样的伪随机数生成器状态。
如果random_state已经是Generator或RandomState实例，则使用该实例。如果random_state是一个整数，则使用一个新的RandomState实例，并使用random_state进行种子化。如果random_state为None（默认），则使用numpy.random.RandomState单例。
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
是否使用‘percentile’自举法（‘percentile’），‘basic’（又名‘reverse’）自举法（‘basic’），或者校正和加速的自举法（‘BCa’，默认值）。
scipy.stats.combine_pvalues
scipy.stats.combine_pvalues(pvalues, method='fisher', weights=None)
结合对同一假设相关的独立测试的 p 值。
这些方法仅用于组合基于连续分布的假设检验的 p 值。
每种方法假设在零假设下，p 值是独立采样且均匀分布于区间 [0, 1]。计算一个检验统计量（每种方法不同），并根据此检验统计量在零假设下的分布计算组合 p 值。
参数：
pvaluesarray_like，1-D
假设来自基于连续分布的独立测试的 p 值数组。
方法
要使用的方法名称来组合 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
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()
scipy.stats.boxcox_normmax
scipy.stats.boxcox_normmax(x, brack=None, method='pearsonr', optimizer=None)
计算输入数据的最佳 Box-Cox 变换参数。
参数：
x array_like
输入数组。所有条目必须为正有限实数。
brack 2-元组，可选，默认为 (-2.0, 2.0)
默认 optimize.brent 求解器进行向下斜坡搜索的起始区间。请注意，这在大多数情况下并不重要；最终结果允许在此区间之外。如果传递了 optimizer，则 brack 必须为 None。
method str，可选
确定最佳变换参数（boxcox 的 lmbda 参数）的方法。选项包括：
‘pearsonr’（默认）
最大化皮尔逊相关系数，使 y = boxcox(x) 与 y 的期望值在 x 若为正态分布时保持一致。
‘mle’
最大化对数似然 boxcox_llf。这是 boxcox 中使用的方法。
‘all’
使用所有可用的优化方法，并返回所有结果。有助于比较不同的方法。
optimizer 可调用函数，可选
optimizer 是一个接受一个参数的可调用函数：
funcallable
要最小化的目标函数。fun 接受一个参数，即 Box-Cox 变换参数 lmbda，并返回在提供的参数处函数值（例如，负对数似然）。optimizer 的任务是找到最小化 fun 的 lmbda 值。
并返回一个对象，如 scipy.optimize.OptimizeResult 的实例，其中保存了属性 x 的最佳 lmbda 值。
更多信息请参见下面的示例或 scipy.optimize.minimize_scalar 的文档。
返回：
maxlog 浮点数或 ndarray
找到的最佳变换参数。对于 method='all'，这是一个数组而不是标量。
参见
boxcox、boxcox_llf、boxcox_normplot、scipy.optimize.minimize_scalar 的示例
示例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
我们可以生成一些数据，并以各种方式确定最佳的 lmbda：
>>> rng = np.random.default_rng()
>>> x = stats.loggamma.rvs(5, size=30, random_state=rng) + 5
>>> y, lmax_mle = stats.boxcox(x)
>>> lmax_pearsonr = stats.boxcox_normmax(x)
>>> lmax_mle
2.217563431465757
>>> lmax_pearsonr
2.238318660200961
>>> stats.boxcox_normmax(x, method='all')
array([2.23831866, 2.21756343])
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> prob = stats.boxcox_normplot(x, -10, 10, plot=ax)
>>> ax.axvline(lmax_mle, color='r')
>>> ax.axvline(lmax_pearsonr, color='g', ls='--')
>>> plt.show()
或者，我们可以定义自己的优化器函数。假设我们只对区间[6, 7]上的lmbda值感兴趣，我们希望使用scipy.optimize.minimize_scalar，设置method='bounded'，并且在优化对数似然函数时希望使用更严格的公差。为此，我们定义一个接受位置参数fun的函数，并使用scipy.optimize.minimize_scalar来最小化fun，同时满足提供的边界和公差：
>>> from scipy import optimize
>>> options = {'xatol': 1e-12} # absolute tolerance on `x`
>>> def optimizer(fun):
... return optimize.minimize_scalar(fun, bounds=(6, 7),
... method="bounded", options=options)
>>> stats.boxcox_normmax(x, optimizer=optimizer)
6.000...
scipy.stats.boxcox_llf
scipy.stats.boxcox_llf(lmb, data)
Box-Cox 对数似然函数。
参数：
lmb 标量
Box-Cox 转换的参数。详情请见 boxcox。
data array_like
计算 Box-Cox 对数似然函数的数据。如果 data 是多维的，则沿着第一轴计算对数似然。
返回：
llf 浮点数或者 ndarray
给定 lmb 的 Box-Cox 对数似然函数。对于 1-D data 是一个浮点数，否则是一个数组。
参见
boxcox, probplot, boxcox_normplot, boxcox_normmax
注意事项
Box-Cox 对数似然函数在这里定义为
[llf = (\lambda - 1) \sum_i(\log(x_i)) - N/2 \log(\sum_i (y_i - \bar{y})² / N),]
其中 y 是经过 Box-Cox 变换的输入数据 x。
示例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
生成一些随机变量并计算它们的 Box-Cox 对数似然值，使用一系列 lmbda 值：
>>> rng = np.random.default_rng()
>>> x = stats.loggamma.rvs(5, loc=10, size=1000, random_state=rng)
>>> lmbdas = np.linspace(-2, 10)
>>> llf = np.zeros(lmbdas.shape, dtype=float)
>>> for ii, lmbda in enumerate(lmbdas):
... llf[ii] = stats.boxcox_llf(lmbda, x)
同时使用 boxcox 找到最优 lmbda 值：
>>> x_most_normal, lmbda_optimal = stats.boxcox(x)
绘制以 lmbda 为函数的对数似然函数图。添加最优 lmbda 作为水平线来确认确实是最优值：
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(lmbdas, llf, 'b.-')
>>> ax.axhline(stats.boxcox_llf(lmbda_optimal, x), color='r')
>>> ax.set_xlabel('lmbda parameter')
>>> ax.set_ylabel('Box-Cox log-likelihood')
现在添加一些概率图，以展示通过 boxcox 变换的数据在对数似然函数最大化时看起来最接近正态分布：
>>> locs = [3, 10, 4] # 'lower left', 'center', 'lower right'
>>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
... xt = stats.boxcox(x, lmbda=lmbda)
... (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
... ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
... ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
... ax_inset.set_xticklabels([])
... ax_inset.set_yticklabels([])
... ax_inset.set_title(r'$\lambda=%1.2f$' % lmbda)
>>> plt.show()
scipy.stats.yeojohnson
scipy.stats.yeojohnson(x, lmbda=None)
返回经 Yeo-Johnson 功率变换后的数据集。
参数：
xndarray
输入数组。应为一维数组。
lmbdafloat，可选
如果 lmbda 为 None，则找到最大化对数似然函数的 lambda，并将其作为第二个输出参数返回。否则，按给定值进行变换。
返回：
yeojohnson：ndarray
经 Yeo-Johnson 功率变换后的数组。
maxlogfloat，可选
如果 lmbda 参数为 None，则第二个返回参数为最大化对数似然函数的 lambda。
参见
probplot，yeojohnson_normplot，yeojohnson_normmax，yeojohnson_llf，boxcox
注：
Yeo-Johnson 变换由以下式给出：
y = ((x + 1)**lmbda - 1) / lmbda, for x >= 0, lmbda != 0
log(x + 1), for x >= 0, lmbda = 0
-((-x + 1)**(2 - lmbda) - 1) / (2 - lmbda), for x < 0, lmbda != 2
-log(-x + 1), for x < 0, lmbda = 2
与 boxcox 不同，yeojohnson 不要求输入数据为正数。
自 1.2.0 版新增。
参考文献
I. Yeo 和 R.A. Johnson，《改善正态性或对称性的新型功率变换家族》，Biometrika 87.4 (2000)：
示例
>>> 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')
我们现在使用 yeojohnson 对数据进行变换，使其最接近正态分布：
>>> ax2 = fig.add_subplot(212)
>>> xt, lmbda = stats.yeojohnson(x)
>>> prob = stats.probplot(xt, dist=stats.norm, plot=ax2)
>>> ax2.set_title('Probplot after Yeo-Johnson transformation')
>>> plt.show()
scipy.stats.yeojohnson_normmax
scipy.stats.yeojohnson_normmax(x, brack=None)
计算最优的 Yeo-Johnson 变换参数。
计算输入数据的最优 Yeo-Johnson 变换参数，使用最大似然估计。
参数：
x 类似数组
输入数组。
brack 2 元组，可选
用于 optimize.brent 的下坡搜索的起始区间。请注意，在大多数情况下这并不关键；最终结果允许超出此区间。如果为 None，则使用 optimize.fminbound 并设置避免溢出的边界。
返回值：
maxlog 浮点数
找到的最优变换参数。
另请参阅
yeojohnson, yeojohnson_llf, yeojohnson_normplot
注意事项
自版本 1.2.0 起新增。
示例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
生成一些数据并确定最优的 lmbda
>>> rng = np.random.default_rng()
>>> x = stats.loggamma.rvs(5, size=30, random_state=rng) + 5
>>> lmax = stats.yeojohnson_normmax(x)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> prob = stats.yeojohnson_normplot(x, -10, 10, plot=ax)
>>> ax.axvline(lmax, color='r')
>>> plt.show()
scipy.stats.yeojohnson_llf
scipy.stats.yeojohnson_llf(lmb, data)
Yeo-Johnson 对数似然函数。
参数：
lmb标量
Yeo-Johnson 变换的参数。详情请参阅yeojohnson。
数据array_like
用于计算 Yeo-Johnson 对数似然的数据。如果data是多维的，则沿第一轴计算对数似然。
返回：
llf浮点数
给定lmb的 Yeo-Johnson 对数似然函数。
另请参阅
yeojohnson，probplot，yeojohnson_normplot，yeojohnson_normmax
注意事项
Yeo-Johnson 对数似然函数在这里定义为
[llf = -N/2 \log(\hat{\sigma}²) + (\lambda - 1) \sum_i \text{ sign }(x_i)\log(|x_i| + 1)]
其中(\hat{\sigma}²)是 Yeo-Johnson 变换后输入数据x的估计方差。
版本 1.2.0 中的新功能。
举例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
生成一些随机变量，并计算它们的 Yeo-Johnson 对数似然值，用一系列lmbda值：
>>> x = stats.loggamma.rvs(5, loc=10, size=1000)
>>> lmbdas = np.linspace(-2, 10)
>>> llf = np.zeros(lmbdas.shape, dtype=float)
>>> for ii, lmbda in enumerate(lmbdas):
... llf[ii] = stats.yeojohnson_llf(lmbda, x)
还可以使用yeojohnson找到最优的lmbda值：
>>> x_most_normal, lmbda_optimal = stats.yeojohnson(x)
绘制对数似然函数作为lmbda的函数。添加最优lmbda作为水平线以检查是否确实是最优：
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(lmbdas, llf, 'b.-')
>>> ax.axhline(stats.yeojohnson_llf(lmbda_optimal, x), color='r')
>>> ax.set_xlabel('lmbda parameter')
>>> ax.set_ylabel('Yeo-Johnson log-likelihood')
现在添加一些概率图，显示对数似然函数最大化的地方，用yeojohnson变换后的数据看起来最接近正态分布：
>>> locs = [3, 10, 4] # 'lower left', 'center', 'lower right'
>>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
... xt = stats.yeojohnson(x, lmbda=lmbda)
... (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
... ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
... ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
... ax_inset.set_xticklabels([])
... ax_inset.set_yticklabels([])
... ax_inset.set_title(r'$\lambda=%1.2f$' % lmbda)
>>> plt.show()
scipy.stats.obrientransform
scipy.stats.obrientransform(*samples)
计算输入数据（任意数量的数组）上的 O’Brien 变换。
用于在运行单因素统计之前测试方差的均匀性。*samples中的每个数组都是因素的一个水平。如果在转换后的数据上运行 f_oneway，并且发现显著性，则方差不相等。来自 Maxwell 和 Delaney [1]，p.112。
参数：
sample1, sample2, …array_like
任意数量的数组。
返回：
obrientransformndarray
用于 ANOVA 的转换数据。结果的第一个维度对应于转换数组的序列。如果给定的数组都是相同长度的 1-D 数组，则返回值是一个 2-D 数组；否则它是一个对象类型的 1-D 数组，其中每个元素都是一个 ndarray。
参考文献
[1]
S. E. Maxwell 和 H. D. Delaney，“Designing Experiments and Analyzing Data: A Model Comparison Perspective”，Wadsworth，1990 年。
示例
我们将测试以下数据集的方差差异。
>>> x = [10, 11, 13, 9, 7, 12, 12, 9, 10]
>>> y = [13, 21, 5, 10, 8, 14, 10, 12, 7, 15]
对数据应用 O’Brien 变换。
>>> from scipy.stats import obrientransform
>>> tx, ty = obrientransform(x, y)
使用 scipy.stats.f_oneway 对转换数据应用单因素 ANOVA 检验。
>>> from scipy.stats import f_oneway
>>> F, p = f_oneway(tx, ty)
>>> p
0.1314139477040335
如果我们要求 p < 0.05 表示显著性，则我们不能断定方差不同。
scipy.stats.sigmaclip
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.sigmaclip.html#scipy.stats.sigmaclip
scipy.stats.sigmaclip(a, low=4.0, high=4.0)
对数组元素执行迭代的 Sigma 切除。
从完整样本开始，移除所有在临界范围之外的元素，即满足以下条件之一的输入数组 c 的所有元素：
c < mean(c) - std(c)*low
c > mean(c) + std(c)*high
迭代继续进行，直到没有元素在（更新后的）范围之外。
参数：
aarray_like
数据数组，如果不是 1-D，则会展平。
lowfloat，可选
Sigma 切除的下限系数。默认为 4。
highfloat，可选
Sigma 切除的上限系数。默认为 4。
返回：
clippedndarray
带有切除元素的输入数组。
lowerfloat
用于切除的下阈值。
upperfloat
用于切除的上阈值。
示例：
>>> import numpy as np
>>> from scipy.stats import sigmaclip
>>> a = np.concatenate((np.linspace(9.5, 10.5, 31),
... np.linspace(0, 20, 5)))
>>> fact = 1.5
>>> c, low, upp = sigmaclip(a, fact, fact)
>>> c
array([ 9.96666667, 10\. , 10.03333333, 10\. ])
>>> c.var(), c.std()
(0.00055555555555555165, 0.023570226039551501)
>>> low, c.mean() - fact*c.std(), c.min()
(9.9646446609406727, 9.9646446609406727, 9.9666666666666668)
>>> upp, c.mean() + fact*c.std(), c.max()
(10.035355339059327, 10.035355339059327, 10.033333333333333)
>>> a = np.concatenate((np.linspace(9.5, 10.5, 11),
... np.linspace(-100, -50, 3)))
>>> c, low, upp = sigmaclip(a, 1.8, 1.8)
>>> (c == np.linspace(9.5, 10.5, 11)).all()
True
scipy.stats.trimboth
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.trimboth.html#scipy.stats.trimboth
scipy.stats.trimboth(a, proportiontocut, axis=0)
从数组的两端切除一部分项目。
从传递的数组的两端切除传递的项目的比例（即，proportiontocut = 0.1，切片左边的 10% 和 右边的 10% 的分数）。修剪的值是最低和最高的值。如果比例导致非整数切片索引，则切片较少（即，保守地切片 proportiontocut）。
参数：
aarray_like
要修剪的数据。
proportiontocutfloat
要修剪的每端的总数据集的比例（范围在 0-1 之间）。
axisint 或 None，可选
数据修剪的轴。默认为 0。如果为 None，则在整个数组 a 上计算。
返回：
outndarray
数组 a 的修剪版本。修剪内容的顺序未定义。
另请参见
trim_mean
示例
创建一个包含 10 个值的数组，并从每端修剪 10% 的值：
>>> import numpy as np
>>> from scipy import stats
>>> a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> stats.trimboth(a, 0.1)
array([1, 3, 2, 4, 5, 6, 7, 8])
输入数组的元素根据值进行修剪，但输出数组未必按顺序排列。
要修剪的比例向下舍入到最接近的整数。例如，从一个包含 10 个值的数组的每端修剪 25% 的值将返回一个包含 6 个值的数组：
>>> b = np.arange(10)
>>> stats.trimboth(b, 1/4).shape
(6,)
可以沿任何轴或整个数组修剪多维数组：
>>> c = [2, 4, 6, 8, 0, 1, 3, 5, 7, 9]
>>> d = np.array([a, b, c])
>>> stats.trimboth(d, 0.4, axis=0).shape
(1, 10)
>>> stats.trimboth(d, 0.4, axis=1).shape
(3, 2)
>>> stats.trimboth(d, 0.4, axis=None).shape
(6,)
scipy.stats.trim1
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.trim1.html#scipy.stats.trim1
scipy.stats.trim1(a, proportiontocut, tail='right', axis=0)
从传递的数组分布的一端切片掉一部分。
如果 proportiontocut = 0.1，则切除分数‘最左侧’或‘最右侧’的 10%得分。修剪较少，如果比例导致非整数切片索引（即保守地切掉 proportiontocut ）。
参数：
a 数组样式
输入数组。
proportiontocut 浮点数
分数是从分布的“左侧”或“右侧”截掉的。
tail {‘left’，‘right’}，可选
默认为‘right’。
axis 整数或无，可选
用于修剪数据的轴。默认为 0。如果为 None，则在整个数组 a 上计算。
返回：
trim1 数组
缩短版本的数组 a。修剪后内容的顺序未定义。
示例
创建包含 10 个值的数组并修剪其最低值的 20%：
>>> import numpy as np
>>> from scipy import stats
>>> a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> stats.trim1(a, 0.2, 'left')
array([2, 4, 3, 5, 6, 7, 8, 9])
请注意，输入数组的元素按值修剪，但输出数组未必排序。
要修剪的比例向下舍入到最接近的整数。例如，从包含 10 个值的数组中修剪 25%的值将返回包含 8 个值的数组：
>>> b = np.arange(10)
>>> stats.trim1(b, 1/4).shape
(8,)
多维数组可以沿任意轴或整个数组进行修剪：
>>> c = [2, 4, 6, 8, 0, 1, 3, 5, 7, 9]
>>> d = np.array([a, b, c])
>>> stats.trim1(d, 0.8, axis=0).shape
(1, 10)
>>> stats.trim1(d, 0.8, axis=1).shape
(3, 2)
>>> stats.trim1(d, 0.8, axis=None).shape
(6,)
scipy.stats.zmap
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.zmap.html#scipy.stats.zmap
scipy.stats.zmap(scores, compare, axis=0, ddof=0, nan_policy='propagate')
计算相对 z-scores。
返回一个 z-score 数组，即标准化为零均值和单位方差的分数，其中均值和方差是从比较数组计算得出的。
参数：
scores：array_like
用于计算 z-scores 的输入。
compare：array_like
用于计算归一化均值和标准差的输入；假设与scores具有相同的维度。
axis：int 或 None，可选
计算compare的均值和方差的轴。默认为 0。如果为 None，则在整个数组scores上计算。
ddof：int，可选
在标准差计算中的自由度校正。默认为 0。
nan_policy：{‘propagate’, ‘raise’, ‘omit’}，可选
定义如何处理compare中 NaN 的出现。‘propagate’返回 NaN，‘raise’引发异常，‘omit’执行计算时忽略 NaN 值。默认为‘propagate’。请注意，当值为‘omit’时，scores中的 NaN 也会传播到输出，但它们不会影响对非 NaN 值计算的 z-scores。
返回：
zscore：array_like
与scores相同形状的 Z-scores。
注意
此函数保留 ndarray 子类，并且还适用于矩阵和掩码数组（它使用asanyarray而不是asarray作为参数）。
示例
>>> from scipy.stats import zmap
>>> a = [0.5, 2.0, 2.5, 3]
>>> b = [0, 1, 2, 3, 4]
>>> zmap(a, b)
array([-1.06066017, 0\. , 0.35355339, 0.70710678])
scipy.stats.zscore
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.zscore.html#scipy.stats.zscore
scipy.stats.zscore(a, axis=0, ddof=0, nan_policy='propagate')
计算 z 分数。
计算样本中每个值相对于样本均值和标准差的 z 分数。
参数：
aarray_like
一个类似数组的对象，包含样本数据。
axisint 或 None，可选
操作的轴。默认为 0。如果为 None，则在整个数组a上计算。
ddofint，可选
在标准差计算中的自由度修正。默认为 0。
nan_policy，可选
定义输入包含 nan 时的处理方式。‘propagate’返回 nan，‘raise’抛出错误，‘omit’执行计算时忽略 nan 值。默认为‘propagate’。注意，当值为‘omit’时，输入中的 nan 也会传播到输出，但它们不会影响计算非 nan 值的 z 分数。
返回：
zscorearray_like
标准化后的 z 分数，按输入数组a的均值和标准差计算。
另请参阅
算术平均
算术标准差
scipy.stats.gzscore
几何标准分数
注释
此函数保留 ndarray 子类，并且还适用于矩阵和掩码数组（它使用asanyarray而不是asarray作为参数）。
参考文献
[1]
“标准分数”，维基百科，zh.wikipedia.org/wiki/%E6%A8%99%E6%BA%96%E5%88%86%E6%95%B8。
[2]
Huck, S. W., Cross, T. L., Clark, S. B，“克服关于 Z 分数的误解”，《教学统计学》，第 8 卷，第 38-40 页，1986 年
示例
>>> import numpy as np
>>> a = np.array([ 0.7972, 0.0767, 0.4383, 0.7866, 0.8091,
... 0.1954, 0.6307, 0.6599, 0.1065, 0.0508])
>>> from scipy import stats
>>> stats.zscore(a)
array([ 1.1273, -1.247 , -0.0552, 1.0923, 1.1664, -0.8559, 0.5786,
0.6748, -1.1488, -1.3324])
沿指定轴计算，使用 n-1 自由度(ddof=1)计算标准差：
>>> b = np.array([[ 0.3148, 0.0478, 0.6243, 0.4608],
... [ 0.7149, 0.0775, 0.6072, 0.9656],
... [ 0.6341, 0.1403, 0.9759, 0.4064],
... [ 0.5918, 0.6948, 0.904 , 0.3721],
... [ 0.0921, 0.2481, 0.1188, 0.1366]])
>>> stats.zscore(b, axis=1, ddof=1)
array([[-0.19264823, -1.28415119, 1.07259584, 0.40420358],
[ 0.33048416, -1.37380874, 0.04251374, 1.00081084],
[ 0.26796377, -1.12598418, 1.23283094, -0.37481053],
[-0.22095197, 0.24468594, 1.19042819, -1.21416216],
[-0.82780366, 1.4457416 , -0.43867764, -0.1792603 ]])
以nan_policy='omit'为例：
>>> x = np.array([[25.11, 30.10, np.nan, 32.02, 43.15],
... [14.95, 16.06, 121.25, 94.35, 29.81]])
>>> stats.zscore(x, axis=1, nan_policy='omit')
array([[-1.13490897, -0.37830299, nan, -0.08718406, 1.60039602],
[-0.91611681, -0.89090508, 1.4983032 , 0.88731639, -0.5785977 ]])
scipy.stats.gzscore
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.gzscore.html#scipy.stats.gzscore
scipy.stats.gzscore(a, *, axis=0, ddof=0, nan_policy='propagate')
计算几何标准分数。
计算样本中每个严格正值的几何 z 分数，相对于几何平均值和标准差。数学上，可以将几何 z 分数评估为：
gzscore = log(a/gmu) / log(gsigma)
其中gmu（或gsigma）是几何平均值（或标准差）。
参数：
a类似数组
样本数据。
axis整数或 None，可选
操作的轴。默认为 0。如果为 None，则在整个数组a上计算。
ddof整数，可选
在标准差计算中的自由度校正。默认为 0。
nan_policy，可选
定义输入包含 NaN 时的处理方式。‘propagate’返回 NaN，‘raise’引发错误，‘omit’执行计算时忽略 NaN 值。默认为‘propagate’。注意，当值为‘omit’时，输入中的 NaN 也会传播到输出，但它们不会影响对非 NaN 值计算的几何 z 分数。
返回：
gzscore类似数组
几何 z 分数，通过输入数组a的几何平均值和几何标准差标准化。
另请参见
gmean
几何平均值
gstd
几何标准差
zscore
标准分数
注意
此函数保留 ndarray 子类，并且也适用于矩阵和掩码数组（它使用asanyarray而不是asarray作为参数）。
新版本中的新功能 1.8。
参考文献
[1]
“几何标准分数”，维基百科，en.wikipedia.org/wiki/Geometric_standard_deviation#Geometric_standard_score。
示例
从对数正态分布中抽样：
>>> import numpy as np
>>> from scipy.stats import zscore, gzscore
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
>>> mu, sigma = 3., 1. # mean and standard deviation
>>> x = rng.lognormal(mu, sigma, size=500)
显示样本的直方图：
>>> fig, ax = plt.subplots()
>>> ax.hist(x, 50)
>>> plt.show()
显示经典 z 分数标准化样本的直方图。分布被重新缩放，但其形状保持不变。
>>> fig, ax = plt.subplots()
>>> ax.hist(zscore(x), 50)
>>> plt.show()
证明几何 z 分数的分布被重新缩放并且准正态：
>>> fig, ax = plt.subplots()
>>> ax.hist(gzscore(x), 50)
>>> plt.show()
scipy.stats.wasserstein_distance
scipy.stats.wasserstein_distance(u_values, v_values, u_weights=None, v_weights=None)
计算两个 1D 分布之间的第一个 Wasserstein 距离。
这个距离也被称为“地球移动距离”，因为它可以看作是将 (u) 转换为 (v) 所需的最小“工作”量，其中“工作”被定义为必须移动的分布权重乘以它必须移动的距离。
新版本 1.0.0 中新增。
参数：
u_values, v_values 数组样式
观察到的值在（经验）分布中。
u_weights, v_weights 数组样式，可选
每个值的权重。如果未指定，每个值将被赋予相同的权重。u_weights（或v_weights）必须与u_values（或v_values）具有相同的长度。如果权重之和与 1 不同，则它仍必须是正的和有限的，以便权重可以归一化为总和为 1。
返回：
distance 浮点数
分布之间计算的距离。
注意事项
分布 (u) 和 (v) 之间的第一个 Wasserstein 距离是：
[l_1 (u, v) = \inf_{\pi \in \Gamma (u, v)} \int_{\mathbb{R} \times \mathbb{R}} |x-y| \mathrm{d} \pi (x, y)]
其中 (\Gamma (u, v)) 是在第一个和第二因子上的边际分布分别为 (u) 和 (v) 的（概率）分布集合。
如果 (U) 和 (V) 是 (u) 和 (v) 的累积分布函数，则此距离也等于：
[l_1(u, v) = \int_{-\infty}^{+\infty} |U-V|]
请参见[2]以证明两个定义的等价性。
输入分布可以是经验的，因此来自样本，其值有效地成为函数的输入，或者它们可以被视为广义函数，此时它们是位于指定值处的 Dirac delta 函数的加权和。
参考文献
[1]
“Wasserstein metric”，en.wikipedia.org/wiki/Wasserstein_metric
[2]
Ramdas, Garcia, Cuturi "On Wasserstein Two Sample Testing and Related Families of Nonparametric Tests" (2015). arXiv:1509.02237.
示例
>>> from scipy.stats import wasserstein_distance
>>> wasserstein_distance([0, 1, 3], [5, 6, 8])
5.0
>>> wasserstein_distance([0, 1], [0, 1], [3, 1], [2, 2])
0.25
>>> wasserstein_distance([3.4, 3.9, 7.5, 7.8], [4.5, 1.4],
... [1.4, 0.9, 3.1, 7.2], [3.2, 3.5])
4.0781331438047861
scipy.stats.energy_distance
scipy.stats.energy_distance(u_values, v_values, u_weights=None, v_weights=None)
计算两个一维分布之间的能量距离。
1.0.0 版新功能。
参数：
u_values, v_valuesarray_like
观察到的（经验）分布中的值。
u_weights, v_weightsarray_like，可选
每个值的权重。如果未指定，则每个值被分配相同的权重。u_weights（v_weights）必须与u_values（v_values）具有相同的长度。如果权重之和不等于 1，则必须仍为正且有限，以便能够将权重归一化为 1。
返回：
distancefloat
计算的分布之间的距离。
注释
两个分布(u)和(v)之间的能量距离，其累积分布函数分别为(U)和(V)，等于：
[D(u, v) = \left( 2\mathbb E|X - Y| - \mathbb E|X - X'| - \mathbb E|Y - Y'| \right)^{1/2}]
其中(X)和(X')（分别(Y)和(Y')）是独立随机变量，其概率分布为(u)（(v)）。
有时，该量的平方被称为“能量距离”（例如在[2]，[4]），但正如[1]和[3]中所指出的那样，仅上述定义符合距离函数（度量）的公理。
如[2]所示，对于一维实值变量，能量距离与 Cramér-von Mises 距离的非分布自由版本相关联：
[D(u, v) = \sqrt{2} l_2(u, v) = \left( 2 \int_{-\infty}^{+\infty} (U-V)² \right)^{1/2}]
注意，普通的 Cramér-von Mises 标准使用距离的无分布版本。详见[2]（第二部分），关于距离两个版本的更多详细信息。
输入分布可以是经验性的，因此来自其值有效作为函数的输入的样本，或者可以视为广义函数，此时它们是位于指定值处的 Dirac δ函数的加权和。
参考文献
[1]
Rizzo, Szekely，“Energy distance.” Wiley Interdisciplinary Reviews: Computational Statistics，8(1):27-38（2015）。
[2] (1,2,3)
Szekely，“E-statistics: The energy of statistical samples.” Bowling Green State University, Department of Mathematics and Statistics, Technical Report 02-16（2002）。
[3]
“Energy distance”，en.wikipedia.org/wiki/Energy_distance
[4]
Bellemare, Danihelka, Dabney, Mohamed, Lakshminarayanan, Hoyer, Munos，“The Cramer Distance as a Solution to Biased Wasserstein Gradients”（2017）。arXiv:1705.10743。
示例
>>> from scipy.stats import energy_distance
>>> energy_distance([0], [2])
2.0000000000000004
>>> energy_distance([0, 8], [0, 8], [3, 1], [2, 2])
1.0000000000000002
>>> energy_distance([0.7, 7.4, 2.4, 6.8], [1.4, 8. ],
... [2.1, 4.2, 7.4, 8. ], [7.6, 8.8])
0.88003340976158217
scipy.stats.rvs_ratio_uniforms
scipy.stats.rvs_ratio_uniforms(pdf, umax, vmin, vmax, size=1, c=0, random_state=None)
使用比例均匀方法从概率密度函数生成随机样本。
自版本 1.12.0 弃用：rvs_ratio_uniforms将在 SciPy 1.15.0 中移除，推荐使用scipy.stats.sampling.RatioUniforms替代。
参数：
pdfcallable
签名为pdf(x)的函数，与分布的概率密度函数成比例。
umax浮点数
u-方向边界矩形的上限。
vmin浮点数
v-方向边界矩形的下限。
vmax浮点数
v-方向边界矩形的上限。
size整数或整数元组，可选
定义随机变量的数量（默认为 1）。
c浮点数，可选。
比例均匀方法的偏移参数，请参见注意事项。默认为 0。
random_state{None, int, numpy.random.Generator,
如果seed为 None（或np.random），则使用numpy.random.RandomState单例。如果seed为 int，则使用新的RandomState实例，并使用seed作为种子。如果seed已经是Generator或RandomState实例，则使用该实例。
返回：
rvsndarray
根据概率密度函数定义的随机变量。
注意事项
请参阅scipy.stats.sampling.RatioUniforms获取文档。
scipy.stats.fit
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.fit.html#scipy.stats.fit
scipy.stats.fit(dist, data, bounds=None, *, guess=None, method='mle', optimizer=<function differential_evolution>)
将离散或连续分布拟合到数据
给定分布、数据和分布参数的边界，返回参数的最大似然估计。
参数：
distscipy.stats.rv_continuous 或 scipy.stats.rv_discrete
表示要适应数据的分布对象。
data1D array_like
要进行分布拟合的数据。如果数据包含任何 np.nan、np.inf 或 -np.inf，拟合方法将引发 ValueError。
boundsdict 或 元组序列，可选
如果是字典，则每个键是分布参数的名称，相应的值是该参数的下界和上界组成的元组。如果分布仅在该参数的有限值范围内定义，那么不需要该参数的条目；例如，某些分布的参数必须在 [0, 1] 区间内。位置 (loc) 和尺度 (scale) 参数的边界是可选的；默认情况下，它们分别固定为 0 和 1。
如果是一个序列，第 i 个元素是分布的第 i 个参数的下界和上界组成的元组。在这种情况下，必须提供所有分布形状参数的边界。可选地，位置和尺度的边界可以跟随分布形状参数。
如果要固定某个形状（例如已知的形状），则下界和上界可以相等。如果用户提供的下界或上界超出了分布定义域的边界，分布域的边界将替换用户提供的值。类似地，必须为整数的参数将被限制为用户提供边界内的整数值。
guessdict 或 array_like，可选
如果是字典，则每个键是分布的参数名称，相应的值是参数值的猜测。
如果是一个序列，第 i 个元素是分布的第 i 个参数的猜测值。在这种情况下，必须提供所有分布形状参数的猜测值。
如果未提供 guess，则不会将决策变量的猜测传递给优化器。如果提供了 guess，则会将任何缺失参数的猜测设置为下界和上界的均值。必须为整数的参数的猜测值将四舍五入为整数值，而位于用户提供边界和分布定义域交集之外的猜测值将被剪裁。
method
使用method="mle"（默认），通过最小化负对数似然函数来计算拟合。对于超出分布支持的观测值，应用大的有限惩罚（而不是无限的负对数似然）。使用method="mse"，通过最小化负对数产品间距函数来计算拟合。对于超出支持的观测值，应用相同的惩罚。我们遵循[1]的方法，该方法适用于具有重复观测样本。
optimizercallable, optional
optimizer是一个可调用对象，接受以下位置参数。
funcallable
要优化的目标函数。fun接受一个参数x，分布的候选形状参数，并返回给定x、dist和提供的data的目标函数值。optimizer的工作是找到最小化fun的决策变量值。
optimizer还必须接受以下关键字参数。
boundssequence of tuples
决策变量值的边界；每个元素将是包含决策变量下限和上限的元组。
如果提供了guess，optimizer还必须接受以下关键字参数。
x0array_like
每个决策变量的猜测。
如果分布有任何必须是整数的形状参数，或者如果分布是离散的且位置参数不固定，optimizer还必须接受以下关键字参数。
integralityarray_like of bools
对于每个决策变量，如果决策变量必须被限制为整数值，则为 True，如果决策变量是连续的，则为 False。
optimizer必须返回一个对象，例如scipy.optimize.OptimizeResult的实例，其中将决策变量的最优值保存在属性x中。如果提供了fun、status或message属性，它们将包含在由fit返回的结果对象中。
返回：
resultFitResult
具有以下字段的对象。
paramsnamedtuple
包含分布形状参数、位置和（如果适用）尺度的最大似然估计的命名元组。
successbool 或 None
优化器是否考虑优化是否成功终止。
messagestr 或 None
优化器提供的任何状态消息。
对象具有以下方法：
nllf(params=None, data=None)
默认情况下，负对数似然函数在给定数据的拟合params处。接受包含分布的替代形状、位置和尺度以及替代数据数组的元组。
plot(ax=None)
将拟合分布的 PDF/PMF 叠加在数据的归一化直方图上。
另请参见
rv_continuous, rv_discrete
注意事项
当用户提供包含最大似然估计的紧密界限时，优化更有可能收敛到最大似然估计。例如，当将二项分布拟合到数据时，每个样本背后的实验数可能已知，在这种情况下，相应的形状参数n可以固定。
参考文献
[1]
邵勇照和 Marjorie G. Hahn. “最大间隔乘积方法：具有强一致性的统一表达。” 伊利诺伊数学期刊 43.3 (1999): 489-499。
例子
假设我们希望将分布拟合到以下数据。
>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> dist = stats.nbinom
>>> shapes = (5, 0.5)
>>> data = dist.rvs(*shapes, size=1000, random_state=rng)
假设我们不知道数据是如何生成的，但我们怀疑它遵循负二项分布，参数为n和p。 (参见scipy.stats.nbinom.) 我们相信参数n小于 30，且参数p必须在区间[0, 1]内。我们将这些信息记录在变量bounds中，并将其传递给fit。
>>> bounds = [(0, 30), (0, 1)]
>>> res = stats.fit(dist, data, bounds)
fit在用户指定的bounds范围内搜索最佳与数据匹配的值（以最大似然估计的意义）。在这种情况下，它找到了与实际生成数据相似的形状值。
>>> res.params
FitParams(n=5.0, p=0.5028157644634368, loc=0.0) # may vary
我们可以通过在数据的归一化直方图上叠加分布的概率质量函数（形状适合数据）来可视化结果。
>>> import matplotlib.pyplot as plt # matplotlib must be installed to plot
>>> res.plot()
>>> plt.show()
注意，n的估计值恰好是整数；这是因为nbinom PMF 的定义域仅包含整数n，而nbinom对象“知道”这一点。nbinom还知道，形状p必须是介于 0 和 1 之间的值。在这种情况下 - 当分布的域对于参数是有限的时候 - 我们不需要为参数指定边界。
>>> bounds = {'n': (0, 30)} # omit parameter p using a `dict`
>>> res2 = stats.fit(dist, data, bounds)
>>> res2.params
FitParams(n=5.0, p=0.5016492009232932, loc=0.0) # may vary
如果我们希望强制分布在n固定为 6 的情况下进行拟合，我们可以将n的上下界都设为 6。然而，请注意，在这种情况下，优化的目标函数值通常会更差（更高）。
>>> bounds = {'n': (6, 6)} # fix parameter `n`
>>> res3 = stats.fit(dist, data, bounds)
>>> res3.params
FitParams(n=6.0, p=0.5486556076755706, loc=0.0) # may vary
>>> res3.nllf() > res.nllf()
True # may vary
注意，前述示例的数值结果是典型的，但可能会有所不同，因为fit使用的默认优化器scipy.optimize.differential_evolution是随机的。然而，我们可以通过定制优化器的设置来确保可复现性 - 或者完全使用不同的优化器 - 使用optimizer参数。
>>> from scipy.optimize import differential_evolution
>>> rng = np.random.default_rng()
>>> def optimizer(fun, bounds, *, integrality):
... return differential_evolution(fun, bounds, strategy='best2bin',
... seed=rng, integrality=integrality)
>>> bounds = [(0, 30), (0, 1)]
>>> res4 = stats.fit(dist, data, bounds, optimizer=optimizer)
>>> res4.params
FitParams(n=5.0, p=0.5015183149259951, loc=0.0)
scipy.stats.ecdf
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ecdf.html#scipy.stats.ecdf
scipy.stats.ecdf(sample)
样本的经验累积分布函数。
经验累积分布函数（ECDF）是样本底层分布的 CDF 的阶梯函数估计。此函数返回表示经验分布函数及其补集经验生存函数的对象。
参数：
sample1D 数组或scipy.stats.CensoredData
除了数组，支持scipy.stats.CensoredData的实例，包括未审查和右审查的观察结果。目前，其他scipy.stats.CensoredData的实例将导致NotImplementedError。
返回：
resECDFResult
具有以下属性的对象。
cdfEmpiricalDistributionFunction
表示样本的经验累积分布函数的对象。
sfEmpiricalDistributionFunction
表示经验生存函数的对象。
cdf和sf属性本身具有以下属性。
quantilesndarray
定义经验 CDF/SF 的样本中的唯一值。
probabilitiesndarray
指数的分位数对应的概率点估计。
和以下方法：
evaluate(x)：
在参数处评估 CDF/SF。
plot(ax)：
在提供的坐标轴上绘制 CDF/SF。
confidence_interval(confidence_level=0.95)：
计算在分位数值周围 CDF/SF 的置信区间。
注意事项
当样本的每个观测是精确测量时，ECDF 在每个观测点[1]处按1/len(sample)递增。
当观测值为下限、上限或上下限时，数据被称为“审查”，sample可以作为scipy.stats.CensoredData的实例提供。
对于右审查数据，ECDF 由 Kaplan-Meier 估计器给出[2]；目前不支持其他形式的审查。
置信区间根据格林伍德公式或更近期的“指数格林伍德”公式计算，如[4]所述。
参考文献
[1] (1,2,3)
Conover, William Jay. 实用非参数统计. Vol. 350. John Wiley & Sons, 1999.
[2]
Kaplan, Edward L., and Paul Meier. “非参数估计来自不完整观测。” 美国统计协会杂志 53.282 (1958): 457-481.
[3]
Goel, Manish Kumar, Pardeep Khanna, and Jugal Kishore. “理解生存分析：Kaplan-Meier 估计。” 国际阿育吠陀研究杂志 1.4 (2010): 274.
[4]
Sawyer, Stanley. “生存分析中的格林伍德和指数格林伍德置信区间。” www.math.wustl.edu/~sawyer/handouts/greenwood.pdf
示例
非截尾数据
如在示例[1]第 79 页中，从一所高中的男生中随机选择了五个男生。他们的一英里跑步时间记录如下。
>>> sample = [6.23, 5.58, 7.06, 6.42, 5.20] # one-mile run times (minutes)
经验分布函数，用来近似样本男生一英里跑步时间的总体分布函数，计算如下。
>>> from scipy import stats
>>> res = stats.ecdf(sample)
>>> res.cdf.quantiles
array([5.2 , 5.58, 6.23, 6.42, 7.06])
>>> res.cdf.probabilities
array([0.2, 0.4, 0.6, 0.8, 1\. ])
要将结果绘制为阶跃函数：
>>> import matplotlib.pyplot as plt
>>> ax = plt.subplot()
>>> res.cdf.plot(ax)
>>> ax.set_xlabel('One-Mile Run Time (minutes)')
>>> ax.set_ylabel('Empirical CDF')
>>> plt.show()
右截尾数据
如在示例[1]第 91 页中，对十个汽车传动带的使用寿命进行了测试。五次测试因测试中的传动带断裂而结束，但其余测试因其他原因结束（例如，研究资金耗尽，但传动带仍然功能良好）。记录了传动带的行驶里程如下。
>>> broken = [77, 47, 81, 56, 80] # in thousands of miles driven
>>> unbroken = [62, 60, 43, 71, 37]
在测试结束时仍然功能良好的传动带的精确寿命时间是未知的，但已知它们超过了记录在unbroken中的值。因此，这些观察被称为“右截尾”，并且使用scipy.stats.CensoredData来表示数据。
>>> sample = stats.CensoredData(uncensored=broken, right=unbroken)
经验生存函数计算如下。
>>> res = stats.ecdf(sample)
>>> res.sf.quantiles
array([37., 43., 47., 56., 60., 62., 71., 77., 80., 81.])
>>> res.sf.probabilities
array([1\. , 1\. , 0.875, 0.75 , 0.75 , 0.75 , 0.75 , 0.5 , 0.25 , 0\. ])
要将结果绘制为阶跃函数：
>>> ax = plt.subplot()
>>> res.cdf.plot(ax)
>>> ax.set_xlabel('Fanbelt Survival Time (thousands of miles)')
>>> ax.set_ylabel('Empirical SF')
>>> plt.show()
scipy.stats.logrank
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.logrank.html#scipy.stats.logrank
scipy.stats.logrank(x, y, alternative='two-sided')
通过 logrank 测试比较两个样本的生存分布。
参数：
x, yarray_like or CensoredData
根据其经验生存函数比较样本。
alternative, optional
定义备择假设。
零假设是两组，如X和Y的生存分布相同。
可用以下备择假设[4]（默认为‘two-sided’）：
-
‘two-sided’: 两组的生存分布不相同。
-
‘less’: 生存组X更受青睐：在某些时间点上，组X的失效率函数小于组Y的失效率函数。
-
‘greater’: survival of group Y is favored: the group X failure rate function is greater than the group Y failure rate function at some times.
返回：
resLogRankResult
一个包含属性的对象：
statisticfloat ndarray
所计算的统计量（如下定义）。其大小是多数其他 logrank 测试实现返回的平方根大小。
pvaluefloat ndarray
测试的计算 p 值。
参见
scipy.stats.ecdf
注意事项
logrank 测试[1]比较观察到的事件数与零假设下预期事件数之间的差异，即两个样本是否从相同分布中抽取。统计量为
[Z_i = \frac{\sum_{j=1}J(O_{i,j}-E_{i,j})}{\sqrt{\sum_{j=1}J V_{i,j}}} \rightarrow \mathcal{N}(0,1)]
where
[E_{i,j} = O_j \frac{N_{i,j}}{N_j}, \qquad V_{i,j} = E_{i,j} \left(\frac{N_j-O_j}{N_j}\right) \left(\frac{N_j-N_{i,j}}{N_j-1}\right),]
(i) denotes the group (i.e. it may assume values (x) or (y), or it may be omitted to refer to the combined sample) (j) denotes the time (at which an event occurred), (N) is the number of subjects at risk just before an event occurred, and (O) is the observed number of events at that time.
由logrank返回的statistic (Z_x)是许多其他实现返回的统计量的（带符号的）平方根。在零假设下，(Z_x**2)渐近地按自由度为一的卡方分布分布。因此，(Z_x)渐近地按标准正态分布分布。使用(Z_x)的优势在于保留了符号信息（即观察到的事件数是否倾向于少于或大于零假设下预期的数量），从而允许scipy.stats.logrank提供单侧备择假设。
参考文献
[1]
Mantel N. “评估生存数据及其相关的两个新秩次统计量。”《癌症化疗报告》，50(3):163-170，PMID: 5910392，1966 年
[2]
Bland, Altman, “对数秩检验”，BMJ，328:1073，DOI:10.1136/bmj.328.7447.1073，2004 年
[3]
“对数秩检验”，维基百科，zh.wikipedia.org/wiki/对数秩检验
[4]
Brown, Mark. “关于对数秩检验方差选择的问题。”《生物统计学》，71.1 (1984): 65-74.
[5]
Klein, John P., 和 Melvin L. Moeschberger.《生存分析：截尾和删节数据的技术》。卷 1230. 纽约：Springer，2003 年。
示例
参考文献[2] 比较了两种不同类型复发性恶性胶质瘤患者的生存时间。下面的样本记录了每位患者参与研究的时间（以周为单位）。由于数据是右截尾的：未截尾的观察对应于观察到的死亡，而截尾的观察对应于患者因其他原因离开研究，因此使用了scipy.stats.CensoredData 类。
>>> from scipy import stats
>>> x = stats.CensoredData(
... uncensored=[6, 13, 21, 30, 37, 38, 49, 50,
... 63, 79, 86, 98, 202, 219],
... right=[31, 47, 80, 82, 82, 149]
... )
>>> y = stats.CensoredData(
... uncensored=[10, 10, 12, 13, 14, 15, 16, 17, 18, 20, 24, 24,
... 25, 28,30, 33, 35, 37, 40, 40, 46, 48, 76, 81,
... 82, 91, 112, 181],
... right=[34, 40, 70]
... )
我们可以计算和可视化两组的经验生存函数如下。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> ax = plt.subplot()
>>> ecdf_x = stats.ecdf(x)
>>> ecdf_x.sf.plot(ax, label='Astrocytoma')
>>> ecdf_y = stats.ecdf(y)
>>> ecdf_x.sf.plot(ax, label='Glioblastoma')
>>> ax.set_xlabel('Time to death (weeks)')
>>> ax.set_ylabel('Empirical SF')
>>> plt.legend()
>>> plt.show()
经验生存函数的视觉检查表明，两组的生存时间倾向于不同。为了正式评估这种差异是否在 1%水平上显著，我们使用了对数秩检验。
>>> res = stats.logrank(x=x, y=y)
>>> res.statistic
-2.73799...
>>> res.pvalue
0.00618...
p 值小于 1%，因此我们可以认为数据证据不支持零假设，支持备择假设，即两个生存函数之间存在差异。
scipy.stats.directional_stats
scipy.stats.directional_stats(samples, *, axis=0, normalize=True)
计算方向数据的样本统计量。
计算样本向量的方向平均值（也称为平均方向向量）和平均结果长度。
方向平均值是向量数据的“首选方向”的度量。它类似于样本均值，但在数据长度无关紧要时使用（例如单位向量）。
平均结果长度是一个介于 0 和 1 之间的值，用于量化方向数据的分散程度：平均结果长度越小，分散程度越大。关于涉及平均结果长度的方向方差的多个定义可见[1]和[2]。
参数：
samplesarray_like
输入数组。必须至少是二维的，且输入的最后一个轴必须与向量空间的维数对应。当输入恰好是二维时，这意味着数据的每一行都是一个向量观测值。
axis整数，默认为 0
计算方向平均值的轴。
normalize: 布尔值，默认为 True
如果为 True，则将输入标准化，以确保每个观测值都是单位向量。如果观测值已经是单位向量，则考虑将其设置为 False，以避免不必要的计算。
返回：
resDirectionalStats
一个包含属性的对象：
mean_directionndarray
方向平均值。
mean_resultant_lengthndarray
平均结果长度 [1]。
另请参阅
circmean
循环均值；即 2D 角度的方向均值。
circvar
循环方差；即 2D 角度的方向方差。
注意事项
此处使用了来自[1]的方向平均值定义。假设观测值是单位向量，则计算如下。
mean = samples.mean(axis=0)
mean_resultant_length = np.linalg.norm(mean)
mean_direction = mean / mean_resultant_length
此定义适用于方向数据（即每个观测的大小无关紧要的向量数据），但不适用于轴向数据（即每个观测的大小和符号都无关紧要的向量数据）。
已经提出了几个涉及平均结果长度 R 的方向方差的定义，包括 1 - R [1]、 1 - R**2 [2] 和 2 * (1 - R) [2]。与选择其中一个不同，此函数返回 R 作为属性 mean_resultant_length，以便用户可以计算其首选的分散度量。
参考
[1] (1,2,3,4)
Mardia, Jupp. (2000). Directional Statistics (p. 163). Wiley.
[2] (1,2,3,4)
en.wikipedia.org/wiki/Directional_statistics
示例
>>> import numpy as np
>>> from scipy.stats import directional_stats
>>> data = np.array([[3, 4], # first observation, 2D vector space
... [6, -8]]) # second observation
>>> dirstats = directional_stats(data)
>>> dirstats.mean_direction
array([1., 0.])
相比之下，向量的常规样本均值会受每个观测值的大小的影响。此外，结果不会是一个单位向量。
>>> data.mean(axis=0)
array([4.5, -2.])
directional_stats的一个典型用例是在球面上寻找一组观测值的有意义中心，例如地理位置。
>>> data = np.array([[0.8660254, 0.5, 0.],
... [0.8660254, -0.5, 0.]])
>>> dirstats = directional_stats(data)
>>> dirstats.mean_direction
array([1., 0., 0.])
另一方面，常规样本均值的结果不位于球面表面上。
>>> data.mean(axis=0)
array([0.8660254, 0., 0.])
该函数还返回平均结果长度，可用于计算方向方差。例如，使用定义 Var(z) = 1 - R 来自于[2]，其中 R 是平均结果长度，我们可以计算上述示例中向量的方向方差为：
>>> 1 - dirstats.mean_resultant_length
0.13397459716167093
scipy.stats.circmean
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.circmean.html#scipy.stats.circmean
scipy.stats.circmean(samples, high=6.283185307179586, low=0, axis=None, nan_policy='propagate', *, keepdims=False)
计算范围内样本的圆均值。
参数：
samples类似数组
输入数组。
highfloat 或 int，可选
样本范围的高边界。默认为 2*pi。
lowfloat 或 int，可选
样本范围的低边界。默认为 0。
axisint 或 None，默认为 None
如果是 int，则为计算统计量的输入轴。输入的每个轴切片（例如行）的统计量将出现在输出的相应元素中。如果为 None，则在计算统计量之前将对输入进行展平。
nan_policy
定义如何处理输入的 NaN。
-
propagate: 如果计算统计量的轴切片（例如行）中存在 NaN，则输出的相应条目将为 NaN。 -
omit: 在执行计算时将忽略 NaN。如果在计算统计量的轴切片中剩余的数据不足，则输出的相应条目将为 NaN。 -
raise: 如果存在 NaN，则会引发ValueError。
keepdimsbool，默认为 False
如果设置为 True，则减少的轴将作为尺寸为一的维度保留在结果中。选择此选项可确保结果正确地广播到输入数组。
返回：
circmeanfloat
圆均值。
另请参阅
circstd
圆标准差。
circvar
圆方差。
注意事项
从 SciPy 1.9 开始，np.matrix 输入（不建议用于新代码）在执行计算之前将转换为 np.ndarray。在这种情况下，输出将是适当形状的标量或 np.ndarray，而不是二维的 np.matrix。同样，虽然忽略掩码数组的屏蔽元素，输出将是标量或 np.ndarray，而不是具有 mask=False 的掩码数组。
示例
为简单起见，所有角度都以度数打印出来。
>>> import numpy as np
>>> from scipy.stats import circmean
>>> import matplotlib.pyplot as plt
>>> angles = np.deg2rad(np.array([20, 30, 330]))
>>> circmean = circmean(angles)
>>> np.rad2deg(circmean)
7.294976657784009
>>> mean = angles.mean()
>>> np.rad2deg(mean)
126.66666666666666
绘制并比较圆均值与算术平均值。
>>> plt.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
... np.sin(np.linspace(0, 2*np.pi, 500)),
... c='k')
>>> plt.scatter(np.cos(angles), np.sin(angles), c='k')
>>> plt.scatter(np.cos(circmean), np.sin(circmean), c='b',
... label='circmean')
>>> plt.scatter(np.cos(mean), np.sin(mean), c='r', label='mean')
>>> plt.legend()
>>> plt.axis('equal')
>>> plt.show()
scipy.stats.circvar
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.circvar.html#scipy.stats.circvar
scipy.stats.circvar(samples, high=6.283185307179586, low=0, axis=None, nan_policy='propagate', *, keepdims=False)
计算假定在范围内的样本的圆形方差。
参数：
samplesarray_like
输入数组。
highfloat 或 int，可选
样本范围的高边界。默认为 2*pi。
lowfloat 或 int，可选
样本范围的低边界。默认为 0。
axisint 或 None，默认为 None
如果是 int，则是计算统计量的输入的轴。输入的每个轴切片（例如行）的统计量将出现在输出的相应元素中。如果是 None，则在计算统计量之前会将输入展平。
nan_policy
定义如何处理输入的 NaN 值。
-
propagate: 如果在计算统计量的轴切片（例如行）中存在 NaN，则输出的相应条目将是 NaN。 -
omit: 在执行计算时将省略 NaN。如果在计算统计量的轴切片中剩余的数据不足，则输出的相应条目将为 NaN。 -
raise: 如果存在 NaN，则会引发ValueError。
keepdimsbool，默认为 False
如果设置为 True，则被缩减的轴会以尺寸为一的维度保留在结果中。使用此选项，结果将正确地与输入数组进行广播。
返回：
circvarfloat
圆形方差。
参见
circmean
圆形平均值。
circstd
圆形标准差。
注意
这里使用的圆形方差的定义是 1-R，其中 R 是平均结果向量。返回的值在范围 [0, 1] 内，0 表示无方差，1 表示大方差。在小角度的极限情况下，该值类似于“线性”方差的一半。
从 SciPy 1.9 开始，np.matrix 输入（不建议新代码使用）在执行计算之前会转换为 np.ndarray。在这种情况下，输出将是一个适当形状的标量或 np.ndarray，而不是一个 2D 的 np.matrix。同样地，忽略掩码数组的掩码元素后，输出将是一个标量或 np.ndarray，而不是带有 mask=False 的掩码数组。
参考文献
[1]
Fisher, N.I. Statistical analysis of circular data. Cambridge University Press, 1993.
示例
>>> import numpy as np
>>> from scipy.stats import circvar
>>> import matplotlib.pyplot as plt
>>> samples_1 = np.array([0.072, -0.158, 0.077, 0.108, 0.286,
... 0.133, -0.473, -0.001, -0.348, 0.131])
>>> samples_2 = np.array([0.111, -0.879, 0.078, 0.733, 0.421,
... 0.104, -0.136, -0.867, 0.012, 0.105])
>>> circvar_1 = circvar(samples_1)
>>> circvar_2 = circvar(samples_2)
绘制样本。
>>> fig, (left, right) = plt.subplots(ncols=2)
>>> for image in (left, right):
... image.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
... np.sin(np.linspace(0, 2*np.pi, 500)),
... c='k')
... image.axis('equal')
... image.axis('off')
>>> left.scatter(np.cos(samples_1), np.sin(samples_1), c='k', s=15)
>>> left.set_title(f"circular variance: {np.round(circvar_1, 2)!r}")
>>> right.scatter(np.cos(samples_2), np.sin(samples_2), c='k', s=15)
>>> right.set_title(f"circular variance: {np.round(circvar_2, 2)!r}")
>>> plt.show()
scipy.stats.circstd
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.circstd.html#scipy.stats.circstd
scipy.stats.circstd(samples, high=6.283185307179586, low=0, axis=None, nan_policy='propagate', *, normalize=False, keepdims=False)
计算假设样本在范围 [low, high] 中的圆形标准偏差。
参数：
samples：array_like
输入数组。
high：float 或 int，可选
样本范围的高边界。默认为 2*pi。
low：float 或 int，可选
样本范围的低边界。默认为 0。
normalize：boolean，可选
如果为 True，则返回的值等于 sqrt(-2*log(R))，并且不依赖于变量单位。如果为 False（默认），返回的值将按 ((high-low)/(2*pi)) 缩放。
axis：int 或 None，默认为 None
如果为整数，则计算统计量的输入轴（例如行）。输入的每个轴切片的统计量将出现在输出的相应元素中。如果为 None，则在计算统计量之前将进行展平。
nan_policy：{‘propagate’, ‘omit’, ‘raise’}
定义如何处理输入的 NaN。
-
propagate：如果在计算统计量的轴切片（例如行）中存在 NaN，则输出的对应条目将为 NaN。 -
omit：执行计算时将忽略 NaN。如果在计算统计量的轴切片中保留的数据不足，则输出的对应条目将为 NaN。 -
raise：如果存在 NaN，则会引发ValueError。
keepdims：bool，默认为 False
如果设置为 True，则减少的轴将作为大小为一的维度保留在结果中。使用此选项，结果将正确广播到输入数组。
返回：
circstd：float
圆形标准偏差。
另请参见
circmean
圆形均值。
circvar
圆形方差。
注意
这使用了来自[1]的圆形标准偏差的定义。本质上，计算如下。
import numpy as np
C = np.cos(samples).mean()
S = np.sin(samples).mean()
R = np.sqrt(C**2 + S**2)
l = 2*np.pi / (high-low)
circstd = np.sqrt(-2*np.log(R)) / l
在小角度极限下，它返回接近‘线性’标准偏差的数字。
从 SciPy 1.9 开始，np.matrix 输入（不推荐新代码使用）在执行计算之前将转换为 np.ndarray。在这种情况下，输出将是一个标量或适当形状的 np.ndarray，而不是 2D 的 np.matrix。类似地，虽然忽略掩码数组的掩码元素，但输出将是一个标量或 np.ndarray，而不是带有 mask=False 的掩码数组。
引用
[1]
Mardia, K. V. (1972). 2. 在 方向数据的统计 (pp. 18-24). Academic Press. DOI:10.1016/C2013-0-07425-7。
示例
>>> import numpy as np
>>> from scipy.stats import circstd
>>> import matplotlib.pyplot as plt
>>> samples_1 = np.array([0.072, -0.158, 0.077, 0.108, 0.286,
... 0.133, -0.473, -0.001, -0.348, 0.131])
>>> samples_2 = np.array([0.111, -0.879, 0.078, 0.733, 0.421,
... 0.104, -0.136, -0.867, 0.012, 0.105])
>>> circstd_1 = circstd(samples_1)
>>> circstd_2 = circstd(samples_2)
绘制样本。
>>> fig, (left, right) = plt.subplots(ncols=2)
>>> for image in (left, right):
... image.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
... np.sin(np.linspace(0, 2*np.pi, 500)),
... c='k')
... image.axis('equal')
... image.axis('off')
>>> left.scatter(np.cos(samples_1), np.sin(samples_1), c='k', s=15)
>>> left.set_title(f"circular std: {np.round(circstd_1, 2)!r}")
>>> right.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
... np.sin(np.linspace(0, 2*np.pi, 500)),
... c='k')
>>> right.scatter(np.cos(samples_2), np.sin(samples_2), c='k', s=15)
>>> right.set_title(f"circular std: {np.round(circstd_2, 2)!r}")
>>> plt.show()
scipy.stats.sobol_indices
scipy.stats.sobol_indices(*, func, n, dists=None, method='saltelli_2010', random_state=None)
Sobol'的全局敏感性指数。
参数：
funccallable 或 dict(str, array_like)
如果func是可调用的，则用于计算 Sobol 指数的函数。其签名必须为：
func(x: ArrayLike) -> ArrayLike
具有形状(d, n)的x和形状(s, n)的输出，其中：
-
d是func的输入维度（输入变量数）， -
s是func的输出维度（输出变量数）， -
n是样本数量（见下文n）。
函数评估值必须是有限的。
如果func是字典，则包含来自三个不同数组的函数评估。键必须为：f_A、f_B和f_AB。f_A和f_B应该具有形状(s, n)，而f_AB应该具有形状(d, s, n)。这是一个高级功能，滥用可能导致分析错误。
nint
用于生成矩阵A和B的样本数。必须是 2 的幂次方。对func进行评估的总点数将为n*(d+2)。
distslist(distributions)，可选
每个参数的分布列表。参数的分布取决于应用程序，并应谨慎选择。假设参数是独立分布的，这意味着它们的值之间没有约束或关系。
分布必须是具有ppf方法的类的实例。
如果func是可调用的，则必须指定，否则将被忽略。
methodCallable 或 str，默认为‘saltelli_2010’
用于计算第一阶和总 Sobol 指数的方法。
如果是可调用的，则其签名必须是：
func(f_A: np.ndarray, f_B: np.ndarray, f_AB: np.ndarray)
-> Tuple[np.ndarray, np.ndarray]
具有形状(s, n)的f_A, f_B和形状(d, s, n)的f_AB。这些数组包含来自三组不同样本的函数评估。输出是形状为(s, d)的第一个和总索引的元组。这是一个高级功能，滥用可能导致分析错误。
random_state{None, int, numpy.random.Generator}, 可选
如果random_state是 int 或 None，则使用np.random.default_rng(random_state)创建一个新的numpy.random.Generator。如果random_state已经是Generator实例，则使用提供的实例。
返回：
resSobolResult
具有属性的对象：
形状为(s, d)的第一阶 Sobol 指数。
第一阶 Sobol 指数。
形状为(s, d)的总阶 Sobol 指数。
总阶 Sobol 指数。
以及方法：
bootstrap(confidence_level: float, n_resamples: int) -> BootstrapSobolResult
一种提供指数置信区间的方法。查看
scipy.stats.bootstrap获取更多详细信息。自助法是在一阶和总阶指数上进行的，并且它们作为属性
first_order和total_order存在于 BootstrapSobolResult 中。
注意
Sobol' 方法 [1], [2] 是一种基于方差的敏感性分析，用于获取每个参数对感兴趣量（QoIs，即 func 的输出）方差的贡献。各自的贡献可以用来排列参数，并通过计算模型的有效（或平均）维度来评估模型的复杂性。
注意
假设参数是独立分布的。每个参数仍然可以遵循任何分布。事实上，分布非常重要，应该与参数的实际分布匹配。
它使用函数方差分解来探索
[\mathbb{V}(Y) = \sum_{i}^{d} \mathbb{V}i (Y) + \sum^{d} \mathbb{V}{ij}(Y) + ... + \mathbb{V}(Y),]
引入条件方差：
[\mathbb{V}i(Y) = \mathbb{\mathbb{V}}[\mathbb{E}(Y|x_i)] \qquad \mathbb{V}(Y) = \mathbb{\mathbb{V}}[\mathbb{E}(Y|x_i x_j)] - \mathbb{V}_i(Y) - \mathbb{V}_j(Y),]
Sobol' 指数表示为
[S_i = \frac{\mathbb{V}i(Y)}{\mathbb{V}[Y]} \qquad S =\frac{\mathbb{V}_{ij}(Y)}{\mathbb{V}[Y]}.]
(S_{i}) 对应于一阶项，评估第 i 个参数的贡献，而 (S_{ij}) 对应于二阶项，说明第 i 和第 j 个参数之间交互的贡献。这些方程可以推广到计算更高阶项；然而，它们的计算成本高昂，并且其解释较为复杂。这就是为什么只提供一阶指数的原因。
总阶指数代表了参数对 QoI 方差的全局贡献，定义如下：
[S_{T_i} = S_i + \sum_j S_{ij} + \sum_{j,k} S_{ijk} + ... = 1 - \frac{\mathbb{V}[\mathbb{E}(Y|x_{\sim i})]}{\mathbb{V}[Y]}.]
一阶指数总和最多为 1，而总阶指数至少为 1。如果没有相互作用，则一阶和总阶指数相等，并且一阶和总阶指数总和为 1。
警告
负的 Sobol' 值是由于数值误差造成的。增加点数 n 应该会有所帮助。
为了进行良好的分析，所需的样本数量随问题的维数增加而增加。例如，对于三维问题，考虑至少 n >= 2**12。模型越复杂，需要的样本就越多。
即使对于纯加法模型，由于数值噪声，指数的总和也可能不为 1。
参考文献
[1]
Sobol, I. M.. “Sensitivity analysis for nonlinear mathematical models.” Mathematical Modeling and Computational Experiment, 1:407-414, 1993.
[2]
Sobol, I. M. (2001). “非线性数学模型的全局敏感性指数及其蒙特卡罗估计。”《数学与计算机仿真》，55(1-3):271-280，DOI:10.1016/S0378-4754(00)00270-6，2001.
[3]
Saltelli, A. “利用模型评估计算敏感性指数的最佳方法。”《计算物理通讯》，145(2):280-297，DOI:10.1016/S0010-4655(02)00280-1，2002.
[4]
Saltelli, A., M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, 和 S. Tarantola. “全局敏感性分析入门。” 2007.
[5]
Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, 和 S. Tarantola. “基于方差的模型输出敏感性分析。总敏感性指数的设计和估计器。”《计算物理通讯》，181(2):259-270，DOI:10.1016/j.cpc.2009.09.018，2010.
[6]
Ishigami, T. 和 T. Homma. “计算机模型不确定性分析中的重要性量化技术。” IEEE，DOI:10.1109/ISUMA.1990.151285，1990。
Examples
以下是对 Ishigami 函数的一个例子 [6]
[Y(\mathbf{x}) = \sin x_1 + 7 \sin² x_2 + 0.1 x_3⁴ \sin x_1,]
其中 (\mathbf{x} \in [-\pi, \pi]³)。该函数表现出强非线性和非单调性。
请注意，Sobol’指数假设样本是独立分布的。在本例中，我们使用每个边缘上的均匀分布。
>>> import numpy as np
>>> from scipy.stats import sobol_indices, uniform
>>> rng = np.random.default_rng()
>>> def f_ishigami(x):
... f_eval = (
... np.sin(x[0])
... + 7 * np.sin(x[1])**2
... + 0.1 * (x[2]**4) * np.sin(x[0])
... )
... return f_eval
>>> indices = sobol_indices(
... func=f_ishigami, n=1024,
... dists=[
... uniform(loc=-np.pi, scale=2*np.pi),
... uniform(loc=-np.pi, scale=2*np.pi),
... uniform(loc=-np.pi, scale=2*np.pi)
... ],
... random_state=rng
... )
>>> indices.first_order
array([0.31637954, 0.43781162, 0.00318825])
>>> indices.total_order
array([0.56122127, 0.44287857, 0.24229595])
可以使用自举法获取置信区间。
>>> boot = indices.bootstrap()
然后，这些信息可以很容易地进行可视化。
>>> import matplotlib.pyplot as plt
>>> fig, axs = plt.subplots(1, 2, figsize=(9, 4))
>>> _ = axs[0].errorbar(
... [1, 2, 3], indices.first_order, fmt='o',
... yerr=[
... indices.first_order - boot.first_order.confidence_interval.low,
... boot.first_order.confidence_interval.high - indices.first_order
... ],
... )
>>> axs[0].set_ylabel("First order Sobol' indices")
>>> axs[0].set_xlabel('Input parameters')
>>> axs[0].set_xticks([1, 2, 3])
>>> _ = axs[1].errorbar(
... [1, 2, 3], indices.total_order, fmt='o',
... yerr=[
... indices.total_order - boot.total_order.confidence_interval.low,
... boot.total_order.confidence_interval.high - indices.total_order
... ],
... )
>>> axs[1].set_ylabel("Total order Sobol' indices")
>>> axs[1].set_xlabel('Input parameters')
>>> axs[1].set_xticks([1, 2, 3])
>>> plt.tight_layout()
>>> plt.show()
Note
默认情况下，scipy.stats.uniform 的支持为 [0, 1]。通过参数 loc 和 scale，可以获得 [loc, loc + scale] 上的均匀分布。
这一结果尤为有趣，因为一阶指数 (S_{x_3} = 0)，而其总体指数为 (S_{T_{x_3}} = 0.244)。这意味着与 (x_3) 的高阶交互作用导致了差异。几乎 25% 的观察方差是由 (x_3) 和 (x_1) 之间的相关性造成的，尽管 (x_3) 本身对 QoI 没有影响。
以下提供了关于该函数的 Sobol’指数的视觉解释。让我们在 ([-\pi, \pi]³) 中生成 1024 个样本，并计算输出的值。
>>> from scipy.stats import qmc
>>> n_dim = 3
>>> p_labels = ['$x_1$', '$x_2$', '$x_3$']
>>> sample = qmc.Sobol(d=n_dim, seed=rng).random(1024)
>>> sample = qmc.scale(
... sample=sample,
... l_bounds=[-np.pi, -np.pi, -np.pi],
... u_bounds=[np.pi, np.pi, np.pi]
... )
>>> output = f_ishigami(sample.T)
现在我们可以根据每个参数绘制输出的散点图。这提供了一种视觉方式来理解每个参数对函数输出的影响。
>>> fig, ax = plt.subplots(1, n_dim, figsize=(12, 4))
>>> for i in range(n_dim):
... xi = sample[:, i]
... ax[i].scatter(xi, output, marker='+')
... ax[i].set_xlabel(p_labels[i])
>>> ax[0].set_ylabel('Y')
>>> plt.tight_layout()
>>> plt.show()
现在 Sobol' 又进了一步：通过给定参数值（黑线），对输出值进行条件计算均值。这对应于术语 (\mathbb{E}(Y|x_i))。对这个术语的方差计算给出 Sobol' 指数的分子。
>>> mini = np.min(output)
>>> maxi = np.max(output)
>>> n_bins = 10
>>> bins = np.linspace(-np.pi, np.pi, num=n_bins, endpoint=False)
>>> dx = bins[1] - bins[0]
>>> fig, ax = plt.subplots(1, n_dim, figsize=(12, 4))
>>> for i in range(n_dim):
... xi = sample[:, i]
... ax[i].scatter(xi, output, marker='+')
... ax[i].set_xlabel(p_labels[i])
... for bin_ in bins:
... idx = np.where((bin_ <= xi) & (xi <= bin_ + dx))
... xi_ = xi[idx]
... y_ = output[idx]
... ave_y_ = np.mean(y_)
... ax[i].plot([bin_ + dx/2] * 2, [mini, maxi], c='k')
... ax[i].scatter(bin_ + dx/2, ave_y_, c='r')
>>> ax[0].set_ylabel('Y')
>>> plt.tight_layout()
>>> plt.show()
看看 (x_3)，均值的方差为零，导致 (S_{x_3} = 0)。但我们可以进一步观察到输出的方差在 (x_3) 的参数值范围内并不是恒定的。这种异方差性可以通过更高阶的交互作用来解释。此外，在 (x_1) 上也能注意到异方差性，这表明 (x_3) 和 (x_1) 之间存在交互作用。在 (x_2) 上，方差似乎是恒定的，因此可以假设与这个参数的交互作用为零。
这种情况在视觉上分析起来相当简单——尽管这只是一种定性分析。然而，当输入参数的数量增加时，这种分析变得不现实，因为很难对高阶项进行结论。因此，使用 Sobol' 指数的好处显而易见。
