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

SciPy 1.12 中文文档(二十五)

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

scipy.stats.pmean

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

scipy.stats.pmean(a, p, *, axis=0, dtype=None, weights=None, nan_policy='propagate', keepdims=False)

沿指定轴计算加权幂均值。

数组 (a_i) 的加权幂均值,相关权重为 (w_i),定义如下:

[\left( \frac{ \sum_{i=1}^n w_i a_i^p }{ \sum_{i=1}^n w_i } \right)^{ 1 / p } , ,]

并且,使用相等的权重,它给出:

[\left( \frac{ 1 }{ n } \sum_{i=1}^n a_i^p \right)^{ 1 / p } , .]

p=0 时,返回几何均值。

这个均值也称为广义均值或 Hölder 均值,不应与 Kolmogorov 广义均值混淆,后者也称为拟算术均值或广义 f-均值 [3]

参数:

a数组样式

输入数组、掩码数组或可转换为数组的对象。

p整数或浮点数

指数。

axis整数或 None,默认为 0

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

dtypedtype,可选

返回数组的类型及累加器,其中元素被求和。若未指定 dtype,则默认为 a 的 dtype,除非 a 的整数 dtype 的精度低于默认平台整数。这种情况下,将使用默认平台整数。

weights数组样式,可选

权重数组可以是 1-D 的(其长度必须是给定 axisa 的大小),或者与 a 的形状相同。默认为 None,即每个值的权重为 1.0。

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdims布尔值,默认为 False

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

返回值:

pmeanndarray,参见上述 dtype 参数。

输出数组,包含幂均值数值。

参见

numpy.average

加权平均

gmean

几何均值

hmean

调和均值

注意

幂均值是在输入数组的单个维度上计算的,默认为 axis=0,或者如果 axis=None,则在数组的所有值上计算。对于整数输入,使用 float64 类型的中间值和返回值。

1.9 版本中的新功能。

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

参考文献

[1]

“广义均值”,维基百科en.wikipedia.org/wiki/Generalized_mean

[2]

Norris, N., “广义均值函数的凸性质”,《数理统计学年刊》,第 8 卷,pp. 118-120,1937 年

[3]

Bullen, P.S., 《均值与它们的不等式手册》,2003 年

示例

>>> from scipy.stats import pmean, hmean, gmean
>>> pmean([1, 4], 1.3)
2.639372938300652
>>> pmean([1, 2, 3, 4, 5, 6, 7], 1.3)
4.157111214492084
>>> pmean([1, 4, 7], -2, weights=[3, 1, 3])
1.4969684896631954 

当 p=-1 时,幂均值等于调和平均数:

>>> pmean([1, 4, 7], -1, weights=[3, 1, 3])
1.9029126213592233
>>> hmean([1, 4, 7], weights=[3, 1, 3])
1.9029126213592233 

当 p=0 时,幂均值定义为几何平均数:

>>> pmean([1, 4, 7], 0, weights=[3, 1, 3])
2.80668351922014
>>> gmean([1, 4, 7], weights=[3, 1, 3])
2.80668351922014 

scipy.stats.kurtosis

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

scipy.stats.kurtosis(a, axis=0, fisher=True, bias=True, nan_policy='propagate', *, keepdims=False)

计算数据集的峰度(Fisher 或 Pearson)。

峰度是四阶中心矩除以方差的平方。如果使用 Fisher 的定义,则从结果中减去 3.0,使正态分布的结果为 0.0。

如果偏差为 False,则使用 k 统计量计算峰度以消除来自有偏矩估计器的偏差。

使用 kurtosistest 查看结果是否接近正态分布。

参数:

a 数组

计算峰度的数据。

axis int 或 None,默认为 0

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

fisher bool,可选

如果为 True,则使用 Fisher 的定义(正态 ==> 0.0)。如果为 False,则使用 Pearson 的定义(正态 ==> 3.0)。

bias bool,可选

如果为 False,则对统计偏差进行校正。

nan_policy {'propagate', 'omit', 'raise'}

定义如何处理输入的 NaN。

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

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

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

keepdims bool,默认为 False

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

返回:

kurtosis 数组

沿轴计算值的峰度,当所有值相等时返回 NaN。

注意

自 SciPy 1.9 开始,np.matrix 输入(不建议新代码使用)在执行计算之前会转换为 np.ndarray。在这种情况下,输出将是适当形状的标量或 np.ndarray,而不是 2D 的 np.matrix。类似地,虽然忽略了掩码数组的掩码元素,但输出将是标量或 np.ndarray,而不是带有 mask=False 的掩码数组。

参考资料

[1]

Zwillinger 和 Kokoska(2000)。《CRC 标准概率和统计表格与公式》。Chapman & Hall:纽约。2000 年。

示例

在 Fisher 的定义中,正态分布的峰度为零。在下面的示例中,峰度接近零,因为它是从数据集而不是连续分布计算得出的。

>>> import numpy as np
>>> from scipy.stats import norm, kurtosis
>>> data = norm.rvs(size=1000, random_state=3)
>>> kurtosis(data)
-0.06928694200380558 

具有较高峰度的分布尾部更重。在费舍尔的定义中,正态分布的峰度值为零,可以作为一个参考点。

>>> import matplotlib.pyplot as plt
>>> import scipy.stats as stats
>>> from scipy.stats import kurtosis 
>>> x = np.linspace(-5, 5, 100)
>>> ax = plt.subplot()
>>> distnames = ['laplace', 'norm', 'uniform'] 
>>> for distname in distnames:
...     if distname == 'uniform':
...         dist = getattr(stats, distname)(loc=-2, scale=4)
...     else:
...         dist = getattr(stats, distname)
...     data = dist.rvs(size=1000)
...     kur = kurtosis(data, fisher=True)
...     y = dist.pdf(x)
...     ax.plot(x, y, label="{}, {}".format(distname, round(kur, 3)))
...     ax.legend() 

拉普拉斯分布的尾部比正态分布更重。均匀分布(具有负峰度)的尾部最细。

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

scipy.stats.mode

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

scipy.stats.mode(a, axis=0, nan_policy='propagate', keepdims=False)

返回传递数组中众数(最常见)值的数组。

如果存在多个这样的值,则仅返回一个。还返回众数箱的计数。

参数:

aarray_like

要查找模式的数字,n 维数组。

axisint 或 None,默认为 0

如果是 int,则是计算输入的轴。输入的每个轴切片(例如行)的统计量将显示在输出的相应元素中。如果为None,则在计算统计量之前将被拉直。

nan_policy

定义如何处理输入 NaN 值。

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

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

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

keepdimsbool,默认为 False

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

返回:

modendarray

众数值数组。

countndarray

每个模式的计数数组。

注释

使用numpy.unique计算众数。在 NumPy 版本 1.21 及之后的版本中,即使具有不同二进制表示的所有 NaN 也被视为等效,并计为同一值的不同实例。

根据惯例,空数组的众数为 NaN,相关计数为零。

从 SciPy 1.9 开始,np.matrix输入(不建议新代码使用)在执行计算之前将被转换为np.ndarray。在这种情况下,输出将是适当形状的标量或np.ndarray,而不是 2D 的np.matrix。同样,尽管被屏蔽数组的屏蔽元素被忽略,输出将是标量或np.ndarray,而不是具有mask=False的屏蔽数组。

示例

>>> import numpy as np
>>> a = np.array([[3, 0, 3, 7],
...               [3, 2, 6, 2],
...               [1, 7, 2, 8],
...               [3, 0, 6, 1],
...               [3, 2, 5, 5]])
>>> from scipy import stats
>>> stats.mode(a, keepdims=True)
ModeResult(mode=array([[3, 0, 6, 1]]), count=array([[4, 2, 2, 1]])) 

要获取整个数组的模式,请指定axis=None

>>> stats.mode(a, axis=None, keepdims=True)
ModeResult(mode=[[3]], count=[[5]])
>>> stats.mode(a, axis=None, keepdims=False)
ModeResult(mode=3, count=5) 

scipy.stats.moment

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

scipy.stats.moment(a, moment=1, axis=0, nan_policy='propagate', *, center=None, keepdims=False)

计算样本的平均值关于均值的第 n 阶矩。

矩是一组点形状的特定定量测量。由于其与偏度和峰度的密切关系,通常用于计算偏度和峰度系数。

参数:

aarray_like

输入数组。

momentint 或 int 的 array_like,可选

返回的中心矩的顺序。默认为 1。

int 或 None,默认:0

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

nan_policy

定义如何处理输入的 NaN。

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

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

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

中心float 或 None,可选

用于计算矩的点。这可以是样本均值、原点或任何其他点。如果None(默认),则计算中心作为样本均值。

keepdimsbool,默认:False

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

返回:

关于center的第 n 阶矩ndarray 或 float

沿给定轴或如果 axis 为 None 则所有值的适当矩。矩计算的分母是观察数,不进行自由度校正。

另请参见

kurtosis, skew, describe

数据样本的第 k 阶矩:

[m_k = \frac{1}{n} \sum_{i = 1}^n (x_i - c)^k]

其中n是样本数,c是计算矩的中心。此函数使用平方的指数计算[1]以提高效率。

请注意,如果a是一个空数组(a.size == 0),则具有一个元素的数组momentmoment.size == 1)将与标量momentnp.isscalar(moment))处理方式相同。这可能会产生意外形状的数组。

自 SciPy 1.9 开始,np.matrix 输入(不建议新代码使用)在执行计算之前会被转换为 np.ndarray。在这种情况下,输出将是适当形状的标量或者 np.ndarray,而不是一个 2D 的 np.matrix。同样,当忽略掩码数组的掩码元素时,输出将是一个标量或者 np.ndarray,而不是带有 mask=False 的掩码数组。

参考文献

[1]

eli.thegreenplace.net/2009/03/21/efficient-integer-exponentiation-algorithms

示例

>>> from scipy.stats import moment
>>> moment([1, 2, 3, 4, 5], moment=1)
0.0
>>> moment([1, 2, 3, 4, 5], moment=2)
2.0 

scipy.stats.expectile

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

scipy.stats.expectile(a, alpha=0.5, *, weights=None)

计算指定水平的期望分位数。

期望分位数在相同方式上是期望的泛化,如分位数是中位数的泛化。水平 alpha = 0.5 处的期望分位数是均值(平均值)。更多细节请参阅注释。

参数:

aarray_like

包含期望分位数的数字的数组。

alphafloat,默认值:0.5

分位数的水平;例如,alpha=0.5 给出了平均值。

weightsarray_like,可选

a 中的值相关联的权重数组。 weights 必须与 a 的形状可广播。默认值为 None,即每个值的权重为 1.0。整数值的权重元素的作用类似于重复相应观察中的 a 那么多次。有关更多详细信息,请参阅注释。

返回值:

expectilendarray

样本的经验分位数在水平 alpha 处。

另请参阅

numpy.mean

算术平均值

numpy.quantile

分位数

注释

通常情况下,具有累积分布函数(CDF)(F) 的随机变量 (X) 的水平 (\alpha) 处的分位数由以下方程的唯一解 (t) 给出:

[\alpha E((X - t)+) = (1 - \alpha) E((t - X)+) ,.]

这里,((x)_+ = \max(0, x)) 是 (x) 的正部分。这个方程也可以等价地写作:

[\alpha \int_t^\infty (x - t)\mathrm{d}F(x) = (1 - \alpha) \int_{-\infty}^t (t - x)\mathrm{d}F(x) ,.]

样本 (a_i)(数组 a)的经验分位数在水平 (\alpha) 处(alpha(数组 weights),它读作 (F_a(x) = \frac{1}{\sum_i w_i} \sum_i w_i 1_{a_i \leq x}),其中指示函数 (1_{A})。这导致了在水平 alpha 处的经验分位数的定义,作为以下方程的唯一解 (t):

[\alpha \sum_{i=1}^n w_i (a_i - t)+ = (1 - \alpha) \sum^n w_i (t - a_i)_+ ,.]

对于 (\alpha=0.5),这简化为加权平均。此外,(\alpha) 越大,分位数的值越大。

最后,水平 (\alpha) 处的期望分位数也可以写成一个最小化问题。通常使用的选择是

[\operatorname{argmin}t E(\lvert 1 - \alpha\rvert(t - X)²) ,.]

参考文献

[1]

W. K. Newey 和 J. L. Powell(1987 年),“非对称最小二乘估计和检验”,《计量经济学》, 55, 819-847。

[2]

T. Gneiting (2009). “Making and Evaluating Point Forecasts,” 美国统计协会杂志, 106, 746 - 762. DOI:10.48550/arXiv.0912.0902

Examples

>>> import numpy as np
>>> from scipy.stats import expectile
>>> a = [1, 4, 2, -1]
>>> expectile(a, alpha=0.5) == np.mean(a)
True
>>> expectile(a, alpha=0.2)
0.42857142857142855
>>> expectile(a, alpha=0.8)
2.5714285714285716
>>> weights = [1, 3, 1, 1] 

scipy.stats.skew

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

scipy.stats.skew(a, axis=0, bias=True, nan_policy='propagate', *, keepdims=False)

计算数据集的样本偏度。

对于正态分布的数据,偏度应该大约为零。对于单峰连续分布,偏度值大于零意味着分布的右尾部分权重更大。函数skewtest可用于确定偏度值是否足够接近零,从统计学角度讲。

参数:

andarray

输入数组。

axis整数或 None,默认值:0

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

bias布尔值,可选

如果为 False,则校正统计偏差。

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdims布尔值,默认值:False

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

返回:

skewnessndarray

沿轴线的值的偏斜度,在所有值相等时返回 NaN。

注意

样本偏斜度被计算为费舍尔-皮尔逊偏斜度系数,即。

[g_1=\frac{m_3}{m_2^{3/2}}]

where

[m_i=\frac{1}{N}\sum_{n=1}N(x[n]-\bar{x})i]

是偏样本(i\texttt{th})中心矩,(\bar{x})是样本均值。如果bias为 False,则校正了偏差并计算出调整后的费舍尔-皮尔逊标准化矩系数,即。

[G_1=\frac{k_3}{k_2^{3/2}}= \frac{\sqrt{N(N-1)}}{N-2}\frac{m_3}{m_2^{3/2}}.]

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

参考文献

[1]

Zwillinger, D. 和 Kokoska, S. (2000). CRC 标准概率和统计表和公式。Chapman & Hall: 纽约。2000 年。第 2.2.24.1 节

示例

>>> from scipy.stats import skew
>>> skew([1, 2, 3, 4, 5])
0.0
>>> skew([2, 8, 0, 4, 1, 9, 9, 0])
0.2650554122698573 

scipy.stats.kstat

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

scipy.stats.kstat(data, n=2, *, axis=None, nan_policy='propagate', keepdims=False)

返回第 n 个 k-统计量(目前 1<=n<=4)。

第 n 个 k-统计量 k_n 是第 n 个累积量(\kappa_n)的唯一对称无偏估计量。

参数:

dataarray_like

输入数组。注意,n 维输入被展平。

n整数,{1, 2, 3, 4},可选

默认值为 2。

axis整数或 None,默认值:None

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

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdims布尔值,默认值:False

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

返回:

kstat浮点数

第 n 个 k-统计量。

另见

kstatvar

返回第 n 个 k-统计量的无偏估计方差

moment

返回样本关于均值的第 n 个中心矩。

注:

对于样本大小 n,前几个 k-统计量为:

[k_{1} = \mu k_{2} = \frac{n}{n-1} m_{2} k_{3} = \frac{ n^{2} } {(n-1) (n-2)} m_{3} k_{4} = \frac{ n^{2} [(n + 1)m_{4} - 3(n - 1) m²_{2}]} {(n-1) (n-2) (n-3)}]

其中(\mu)是样本均值,(m_2)是样本方差,(m_i)是第 i 个样本中心矩。

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

参考文献

mathworld.wolfram.com/k-Statistic.html

mathworld.wolfram.com/Cumulant.html

示例

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

随着样本大小的增加,第 n 个矩和第 n 个 k-统计量收敛到相同的数值(尽管它们不完全相同)。在正态分布的情况下,它们收敛到零。

>>> for n in [2, 3, 4, 5, 6, 7]:
...     x = rng.normal(size=10**n)
...     m, k = stats.moment(x, 3), stats.kstat(x, 3)
...     print("%.3g  %.3g  %.3g" % (m, k, m-k))
-0.631 -0.651 0.0194  # random
0.0282 0.0283 -8.49e-05
-0.0454 -0.0454 1.36e-05
7.53e-05 7.53e-05 -2.26e-09
0.00166 0.00166 -4.99e-09
-2.88e-06 -2.88e-06 8.63e-13 

scipy.stats.kstatvar

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

scipy.stats.kstatvar(data, n=2, *, axis=None, nan_policy='propagate', keepdims=False)

返回 k-统计量方差的无偏估计器。

查看kstat以获取 k-统计量的更多详细信息。

参数:

dataarray_like

输入数组。请注意,n 维输入会被展平。

nint,{1, 2},可选

默认为 2。

axisint 或 None,默认值:None

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

nan_policy

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

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

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

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

keepdimsbool,默认值:False

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

返回:

kstatvarfloat

第 n 个 k-统计量的方差。

另请参阅

kstat

返回第 n 个 k-统计量。

moment

返回样本关于均值的第 n 个中心矩。

注意事项

前几个 k-统计量的方差为:

[var(k_{1}) = \frac{\kappa²}{n} var(k_{2}) = \frac{\kappa⁴}{n} + \frac{2\kappa²_{2}}{n - 1} var(k_{3}) = \frac{\kappa⁶}{n} + \frac{9 \kappa_2 \kappa_4}{n - 1} + \frac{9 \kappa²_{3}}{n - 1} + \frac{6 n \kappa³_{2}}{(n-1) (n-2)} var(k_{4}) = \frac{\kappa⁸}{n} + \frac{16 \kappa_2 \kappa_6}{n - 1} + \frac{48 \kappa_{3} \kappa_5}{n - 1} + \frac{34 \kappa²_{4}}{n-1} + \frac{72 n \kappa²_{2} \kappa_4}{(n - 1) (n - 2)} + \frac{144 n \kappa_{2} \kappa²_{3}}{(n - 1) (n - 2)} + \frac{24 (n + 1) n \kappa⁴_{2}}{(n - 1) (n - 2) (n - 3)}]

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

scipy.stats.tmean

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

scipy.stats.tmean(a, limits=None, inclusive=(True, True), axis=None, *, nan_policy='propagate', keepdims=False)

计算修剪均值。

此函数找到给定值的算术平均值,忽略limits外的值。

参数:

a类似数组

数组的值。

limitsNone 或(下限,上限),可选

输入数组中小于下限或大于上限的值将被忽略。当 limits 为 None(默认值)时,使用所有值。元组中的任一限值也可以是 None,表示半开区间。

inclusive(布尔值,布尔值),可选

元组包含(下限标志,上限标志)。这些标志确定是否包括等于下限或上限的值。默认值为(True,True)。

axis整数或 None,默认为:None

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

nan_policy

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

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

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

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

keepdims布尔值,默认为:False

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

返回:

tmeanndarray

修剪均值。

另请参见

trim_mean

返回修剪了两侧比例后的均值。

注释

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

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.tmean(x)
9.5
>>> stats.tmean(x, (3,17))
10.0 

scipy.stats.tvar

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

scipy.stats.tvar(a, limits=None, inclusive=(True, True), axis=0, ddof=1, *, nan_policy='propagate', keepdims=False)

计算修剪的方差。

此函数计算值数组的样本方差,同时忽略超出给定限制的值。

参数:

aarray_like

值数组。

limitsNone 或(下限, 上限),可选

输入数组中小于下限或大于上限的值将被忽略。当 limits 为 None 时,所有值都被使用。元组中的任一限制值也可以为 None,表示半开区间。默认值为 None。

inclusive(bool, bool),可选

一个由(下限标志,上限标志)组成的元组。这些标志确定是否包括等于下限或上限的值。默认值为(True, True)。

axisint 或 None,默认值:0

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

ddofint,可选

自由度的增量。默认值为 1。

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdimsbool,默认值:False

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

返回:

tvarfloat

修剪方差。

注意

tvar计算无偏样本方差,即使用修正因子n / (n - 1)

从 SciPy 1.9 开始,np.matrix输入(不建议在新代码中使用)在执行计算之前将被转换为np.ndarray。在这种情况下,输出将是适当形状的标量或np.ndarray,而不是 2D 的np.matrix。类似地,忽略掩码数组的掩码元素,输出将是适当形状的标量或np.ndarray,而不是带有mask=False的掩码数组。

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.tvar(x)
35.0
>>> stats.tvar(x, (3,17))
20.0 

scipy.stats.tmin

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

scipy.stats.tmin(a, lowerlimit=None, axis=0, inclusive=True, nan_policy='propagate', *, keepdims=False)

计算修剪后的最小值。

此函数沿指定轴找到数组a的最小值,但仅考虑大于指定下限的值。

参数:

aarray_like

值数组。

lowerlimitNone 或浮点数,可选

输入数组中小于给定限制的值将被忽略。当 lowerlimit 为 None 时,将使用所有值。默认值为 None。

axis整数或 None,默认值:0

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

inclusive,可选

此标志确定是否包括与下限完全相等的值。默认值为 True。

nan_policy

定义如何处理输入 NaN。

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

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

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

keepdimsbool,默认值:False

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

返回:

tmin浮点数、整数或 ndarray

修剪后的最小值。

注意事项

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

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.tmin(x)
0 
>>> stats.tmin(x, 13)
13 
>>> stats.tmin(x, 13, inclusive=False)
14 

scipy.stats.tmax

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

scipy.stats.tmax(a, upperlimit=None, axis=0, inclusive=True, nan_policy='propagate', *, keepdims=False)

计算被修剪的最大值。

此函数计算沿给定轴的数组的最大值,同时忽略大于指定上限的值。

参数:

aarray_like

值的数组。

upperlimitNone 或 float,可选

输入数组中大于给定限制的值将被忽略。当 upperlimit 为 None 时,将使用所有值。默认值为 None。

axisint 或 None,默认:0

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

inclusive,可选

此标志确定是否包括等于上限的值。默认值为 True。

nan_policy

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

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

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

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

keepdimsbool,默认值:False

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

返回:

tmaxfloat、int 或 ndarray

被修剪的最大值。

注:

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

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.tmax(x)
19 
>>> stats.tmax(x, 13)
13 
>>> stats.tmax(x, 13, inclusive=False)
12 

scipy.stats.tstd

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

scipy.stats.tstd(a, limits=None, inclusive=(True, True), axis=0, ddof=1, *, nan_policy='propagate', keepdims=False)

计算修剪样本标准差。

此函数找到给定值的样本标准差,忽略给定 limits 外的值。

参数:

a array_like

值数组。

limits None 或(下限,上限),可选

输入数组中小于下限或大于上限的值将被忽略。当限制为 None 时,所有值都被使用。元组中的任一限制值也可以为 None,表示半开区间。默认值为 None。

inclusive(布尔值,布尔值),可选

由(较低标志,较高标志)组成的元组。这些标志确定是否包含值等于下限或上限。默认值为(True,True)。

axis整数或 None,默认值:0

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

ddof整数,可选

自由度的 Delta。默认为 1。

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdims 布尔值,默认值:False

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

返回:

tstd 浮点数

修剪样本标准差。

注意

tstd计算无偏样本标准差,即使用校正因子 n / (n - 1)

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

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.tstd(x)
5.9160797830996161
>>> stats.tstd(x, (3,17))
4.4721359549995796 

scipy.stats.tsem

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

scipy.stats.tsem(a, limits=None, inclusive=(True, True), axis=0, ddof=1, *, nan_policy='propagate', keepdims=False)

计算剪裁后的平均标准误差。

此函数找到给定值的平均标准误差,忽略超出给定limits的值。

参数:

aarray_like

值数组。

limitsNone 或者 (下限, 上限),可选项

输入数组中小于下限或大于上限的值将被忽略。当 limits 为 None 时,将使用所有值。元组中的任何一个限制值也可以是 None,表示半开区间。默认值为 None。

inclusive(布尔值,布尔值),可选项

一个元组,包含(下限标志,上限标志)。这些标志确定是否包括与下限或上限完全相等的值。默认值为(True,True)。

axis整数或 None,默认为 0

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

ddof整数,可选项

自由度增量。默认值为 1。

nan_policy

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

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

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

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

keepdims布尔值,默认为 False

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

返回:

tsemfloat

剪裁后的平均标准误差。

注释

tsem 使用无偏样本标准差,即使用校正因子 n / (n - 1)

从 SciPy 1.9 开始,np.matrix 输入(不建议用于新代码)在执行计算之前将转换为 np.ndarray。在这种情况下,输出将是适当形状的标量或 np.ndarray,而不是 2D 的 np.matrix。类似地,忽略掩码数组的屏蔽元素,输出将是标量或 np.ndarray,而不是具有 mask=False 的掩码数组。

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.tsem(x)
1.3228756555322954
>>> stats.tsem(x, (3,17))
1.1547005383792515 

scipy.stats.variation

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

scipy.stats.variation(a, axis=0, nan_policy='propagate', ddof=0, *, keepdims=False)

计算变异系数。

变异系数是标准偏差除以均值。此函数等效于:

np.std(x, axis=axis, ddof=ddof) / np.mean(x) 

ddof的默认值为 0,但是许多变异系数的定义使用样本标准偏差的无偏样本方差的平方根,对应于ddof=1

函数不取数据均值的绝对值,因此如果均值为负,则返回值为负。

参数:

aarray_like

输入数组。

axisint 或 None,默认值:0

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

nan_policy

定义如何处理输入的 NaN。

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

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

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

ddofint,可选

提供了在计算标准偏差时使用的“Delta Degrees Of Freedom”(自由度)。在计算标准偏差时使用的除数是N - ddof,其中N是元素的数量。ddof必须小于N;如果不是,则结果将是naninf,这取决于N和数组中的值。默认情况下,ddof为零以确保向后兼容性,但建议使用ddof=1以确保计算样本标准偏差作为无偏样本方差的平方根。

keepdimsbool,默认值:False

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

返回:

variationndarray

请求轴上计算的计算变异。

注意事项

处理多种边缘情况而不生成警告:

  • 如果均值和标准偏差都为零,则返回nan

  • 如果均值为零且标准偏差不为零,则返回inf

  • 如果输入长度为零(因为数组长度为零,或所有输入值都是nannan_policy'omit'),则返回nan

  • 如果输入包含inf,则返回nan

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

参考资料

[1]

Zwillinger, D. 和 Kokoska, S.(2000)。CRC 标准概率和统计表格与公式。Chapman & Hall:纽约。2000 年。

示例

>>> import numpy as np
>>> from scipy.stats import variation
>>> variation([1, 2, 3, 4, 5], ddof=1)
0.5270462766947299 

计算包含少量nan值的数组沿给定维度的变化:

>>> x = np.array([[  10.0, np.nan, 11.0, 19.0, 23.0, 29.0, 98.0],
...               [  29.0,   30.0, 32.0, 33.0, 35.0, 56.0, 57.0],
...               [np.nan, np.nan, 12.0, 13.0, 16.0, 16.0, 17.0]])
>>> variation(x, axis=1, ddof=1, nan_policy='omit')
array([1.05109361, 0.31428986, 0.146483  ]) 

scipy.stats.find_repeats

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

scipy.stats.find_repeats(arr)

查找重复项和重复计数。

参数:

arrarray_like

输入数组。此数组被转换为 float64 类型。

返回:

valuesndarray

来自(扁平化的)输入的唯一值,它们是重复的。

countsndarray

相应的“value”重复的次数。

笔记

在 numpy >= 1.9 中,numpy.unique 提供类似的功能。主要区别在于 find_repeats 只返回重复的值。

示例

>>> from scipy import stats
>>> stats.find_repeats([2, 1, 2, 3, 2, 2, 5])
RepeatedResults(values=array([2.]), counts=array([4])) 
>>> stats.find_repeats([[10, 20, 1, 2], [5, 5, 4, 4]])
RepeatedResults(values=array([4.,  5.]), counts=array([2, 2])) 

scipy.stats.rankdata

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

scipy.stats.rankdata(a, method='average', *, axis=None, nan_policy='propagate')

分配排名给数据,适当处理并列值。

默认情况下(axis=None),数据数组首先被展平,返回一个平坦的排名数组。如果需要,可以单独将排名数组重塑为数据数组的形状(请参见示例)。

排名从 1 开始。method参数控制如何对等值分配排名。详细讨论排名方法,请参见[1]

参数:

aarray_like

要排名的值数组。

method,可选

用于对并列元素分配排名的方法。提供以下方法(默认为‘average’):

  • ‘average’:将所有并列值分配的排名的平均值分配给每个值。
  • ‘min’:将所有并列值分配的排名的最小值分配给每个值。(这也称为“竞争”排名。)
  • ‘max’:将所有并列值分配的排名的最大值分配给每个值。
  • ‘dense’:类似于‘min’,但是将下一个最高元素的排名分配给紧接在并列元素之后的排名。
  • ‘ordinal’:所有值都被赋予不同的排名,对应于它们在a中出现的顺序。

axis,可选

执行排名的轴。如果为None,则首先展平数据数组。

nan_policy,可选

定义输入包含 nan 时的处理方式。提供以下选项(默认为‘propagate’):

  • ‘propagate’:通过排名计算传播 nan
  • ‘omit’:在执行排名时忽略 nan 值
  • ‘raise’:引发错误

注意

nan_policy为‘propagate’时,输出是所有 nan 的数组,因为输入中 nan 的排名是未定义的。当nan_policy为‘omit’时,排名其他值时会忽略a中的 nan,并且输出的对应位置是 nan。

版本 1.10 中的新增功能。

返回:

ranksndarray

一个大小与a相同的数组,包含排名分数。

参考文献

[1]

“排名”,en.wikipedia.org/wiki/Ranking

示例

>>> import numpy as np
>>> from scipy.stats import rankdata
>>> rankdata([0, 2, 3, 2])
array([ 1\. ,  2.5,  4\. ,  2.5])
>>> rankdata([0, 2, 3, 2], method='min')
array([ 1,  2,  4,  2])
>>> rankdata([0, 2, 3, 2], method='max')
array([ 1,  3,  4,  3])
>>> rankdata([0, 2, 3, 2], method='dense')
array([ 1,  2,  3,  2])
>>> rankdata([0, 2, 3, 2], method='ordinal')
array([ 1,  2,  4,  3])
>>> rankdata([[0, 2], [3, 2]]).reshape(2,2)
array([[1\. , 2.5],
 [4\. , 2.5]])
>>> rankdata([[0, 2, 2], [3, 2, 5]], axis=1)
array([[1\. , 2.5, 2.5],
 [2\. , 1\. , 3\. ]])
>>> rankdata([0, 2, 3, np.nan, -2, np.nan], nan_policy="propagate")
array([nan, nan, nan, nan, nan, nan])
>>> rankdata([0, 2, 3, np.nan, -2, np.nan], nan_policy="omit")
array([ 2.,  3.,  4., nan,  1., nan]) 

scipy.stats.tiecorrect

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

scipy.stats.tiecorrect(rankvals)

Mann-Whitney U 和 Kruskal-Wallis H 检验的校正系数。

参数:

rankvalsarray_like

一个一维排名序列。通常这将是由rankdata返回的数组。

返回:

factorfloat

U 或 H 的校正因子。

另请参阅

rankdata

为数据分配排名

mannwhitneyu

Mann-Whitney 秩和检验

kruskal

Kruskal-Wallis H 检验

参考文献

[1]

Siegel, S. (1956)《行为科学的非参数统计》。纽约:麦格劳-希尔。

示例

>>> from scipy.stats import tiecorrect, rankdata
>>> tiecorrect([1, 2.5, 2.5, 4])
0.9
>>> ranks = rankdata([1, 3, 2, 4, 5, 7, 2, 8, 4])
>>> ranks
array([ 1\. ,  4\. ,  2.5,  5.5,  7\. ,  8\. ,  2.5,  9\. ,  5.5])
>>> tiecorrect(ranks)
0.9833333333333333 

scipy.stats.trim_mean

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

scipy.stats.trim_mean(a, proportiontocut, axis=0)

返回修剪了分布两端后数组的均值。

如果proportiontocut = 0.1,则切掉分数的‘最左端’和‘最右端’各 10%。切片前对输入进行排序。如果比例导致非整数切片索引,则保守地切掉proportiontocut

参数:

aarray_like

输入数组。

proportiontocutfloat

分布两端要切掉的分数比例。

axisint 或 None,可选

计算修剪均值的轴。默认为 0。如果为 None,则在整个数组a上计算。

返回:

trim_meanndarray

修剪后数组的均值。

参见

trimboth

tmean

计算在给定limits外忽略的修剪均值。

示例

>>> import numpy as np
>>> from scipy import stats
>>> x = np.arange(20)
>>> stats.trim_mean(x, 0.1)
9.5
>>> x2 = x.reshape(5, 4)
>>> x2
array([[ 0,  1,  2,  3],
 [ 4,  5,  6,  7],
 [ 8,  9, 10, 11],
 [12, 13, 14, 15],
 [16, 17, 18, 19]])
>>> stats.trim_mean(x2, 0.25)
array([  8.,   9.,  10.,  11.])
>>> stats.trim_mean(x2, 0.25, axis=1)
array([  1.5,   5.5,   9.5,  13.5,  17.5]) 

scipy.stats.gstd

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

scipy.stats.gstd(a, axis=0, ddof=1)

计算数组的几何标准偏差。

几何标准偏差描述了首选几何平均值的一组数字的扩展。它是一个乘法因子,因此是一个无量纲的量。

定义为log(a)的标准偏差的指数。数学上,人口几何标准偏差可以计算为:

gstd = exp(std(log(a))) 

新版本 1.3.0 中。

参数:

aarray_like

一个类似数组的对象,包含样本数据。

axisint、元组或无,可选

沿其操作的轴。默认为 0。如果为 None,则在整个数组a上计算。

ddofint,可选

在计算几何标准偏差时需要使用自由度修正。默认值为 1。

返回:

gstdndarray 或浮点数

一个几何标准偏差的数组。如果axis为 None 或a是 1 维数组,则返回一个浮点数。

参见

gmean

几何平均数

numpy.std

标准偏差

gzscore

几何标准分数

注释

由于计算需要使用对数,几何标准偏差仅支持严格正值。任何非正或无限值都会引发ValueError。几何标准偏差有时会与标准偏差的指数exp(std(a))混淆。实际上,几何标准偏差是exp(std(log(a)))ddof的默认值与其他包含 ddof 函数的默认值(0)不同,如np.stdnp.nanstd

参考文献

[1]

“几何标准偏差”,维基百科en.wikipedia.org/wiki/Geometric_standard_deviation.

[2]

Kirkwood,T.B.,“几何平均数和离散度度量”,生物统计学,第 35 卷,第 908-909 页,1979 年

示例

找到对数正态分布样本的几何标准偏差。注意,分布的标准偏差为 1,在对数尺度上大约为exp(1)

>>> import numpy as np
>>> from scipy.stats import gstd
>>> rng = np.random.default_rng()
>>> sample = rng.lognormal(mean=0, sigma=1, size=1000)
>>> gstd(sample)
2.810010162475324 

计算多维数组和给定轴的几何标准偏差。

>>> a = np.arange(1, 25).reshape(2, 3, 4)
>>> gstd(a, axis=None)
2.2944076136018947
>>> gstd(a, axis=2)
array([[1.82424757, 1.22436866, 1.13183117],
 [1.09348306, 1.07244798, 1.05914985]])
>>> gstd(a, axis=(1,2))
array([2.12939215, 1.22120169]) 

几何标准偏差进一步处理了掩码数组。

>>> a = np.arange(1, 25).reshape(2, 3, 4)
>>> ma = np.ma.masked_where(a > 16, a)
>>> ma
masked_array(
 data=[[[1, 2, 3, 4],
 [5, 6, 7, 8],
 [9, 10, 11, 12]],
 [[13, 14, 15, 16],
 [--, --, --, --],
 [--, --, --, --]]],
 mask=[[[False, False, False, False],
 [False, False, False, False],
 [False, False, False, False]],
 [[False, False, False, False],
 [ True,  True,  True,  True],
 [ True,  True,  True,  True]]],
 fill_value=999999)
>>> gstd(ma, axis=2)
masked_array(
 data=[[1.8242475707663655, 1.2243686572447428, 1.1318311657788478],
 [1.0934830582350938, --, --]],
 mask=[[False, False, False],
 [False,  True,  True]],
 fill_value=999999) 

scipy.stats.iqr

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

scipy.stats.iqr(x, axis=None, rng=(25, 75), scale=1.0, nan_policy='propagate', interpolation='linear', keepdims=False)

计算沿指定轴的数据的四分位距。

四分位距(IQR)是数据的第 75 百分位数和第 25 百分位数之间的差异。它是一种类似于标准差或方差的离散度量,但对异常值更为稳健 [2]

参数rng允许此函数计算除了实际 IQR 之外的其他百分位范围。例如,设置rng=(0, 100)等效于numpy.ptp

空数组的 IQR 为 np.nan

从版本 0.18.0 开始。

参数:

xarray_like

输入数组或可转换为数组的对象。

axisint 或 None,默认值:None

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

rng两个浮点数的序列,范围在[0,100]之间,可选

要计算范围的百分位数。每个必须在 0 到 100 之间,包括 0 和 100。默认为真实 IQR:(25, 75)。元素的顺序不重要。

scale标量或字符串或实数组成的 array_like,可选

scale 的数值将除以最终结果。也识别以下字符串值:

  • ‘normal’:按(2 \sqrt{2} erf^{-1}(\frac{1}{2}) \approx 1.349)缩放。

默认为 1.0。也允许具有实数 dtype 的 array-like scale,只要它正确广播到输出,使得out / scale是有效的操作。输出的维度取决于输入数组 xaxis 参数和 keepdims 标志。

nan_policy

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

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

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

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

interpolation字符串,可选

指定在百分位边界位于两个数据点ij之间时要使用的插值方法。可用以下选项(默认为‘linear’):

  • ‘linear’: i + (j - i)*fraction,其中fraction是由ij包围的索引的分数部分。
  • ‘lower’: i.
  • ‘higher’: j.
  • ‘nearest’: ij中最近的一个。
  • ‘midpoint’: (i + j)/2.

对于 NumPy >= 1.22.0,numpy.percentilemethod 关键字提供的附加选项也是有效的。

keepdims 布尔值,默认值:False

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

返回:

iqr 标量或 ndarray

如果 axis=None,则返回标量。如果输入包含小于 np.float64 的整数或浮点数,则输出数据类型为 np.float64。否则,输出数据类型与输入相同。

另请参阅

numpy.std, numpy.var

注意事项

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

参考资料

[1]

“四分位距” zh.wikipedia.org/wiki/%E5%9B%9B%E5%88%86%E4%BD%8D%E8%B7%9D

[2]

“尺度的稳健测度” zh.wikipedia.org/wiki/%E5%B0%BA%E5%BA%A6%E7%9A%84%E7%A8%B3%E5%81%A5%E6%B5%8B%E5%BA%A6

[3]

“分位数” zh.wikipedia.org/wiki/%E5%88%86%E4%BD%8D%E6%95%B0

示例

>>> import numpy as np
>>> from scipy.stats import iqr
>>> x = np.array([[10, 7, 4], [3, 2, 1]])
>>> x
array([[10,  7,  4],
 [ 3,  2,  1]])
>>> iqr(x)
4.0
>>> iqr(x, axis=0)
array([ 3.5,  2.5,  1.5])
>>> iqr(x, axis=1)
array([ 3.,  1.])
>>> iqr(x, axis=1, keepdims=True)
array([[ 3.],
 [ 1.]]) 

scipy.stats.sem

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

scipy.stats.sem(a, axis=0, ddof=1, nan_policy='propagate', *, keepdims=False)

计算均值的标准误差。

计算输入数组中值的均值标准误差(或测量标准误差)。

参数:

aarray_like

包含标准误差值的数组。

axisint 或 None,默认值:0

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

ddofint, optional

Delta 自由度。在有限样本中相对于总体方差估计进行偏差调整的自由度数量。默认为 1。

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdimsbool,默认值:False

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

返回:

sndarray 或 float

样本中的均值标准误差,沿着输入轴。

注释

ddof 的默认值与其他包含 ddof 的例程(例如 np.std 和 np.nanstd)使用的默认值(0)不同。

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

示例

沿第一个轴找到标准误差:

>>> import numpy as np
>>> from scipy import stats
>>> a = np.arange(20).reshape(5,4)
>>> stats.sem(a)
array([ 2.8284,  2.8284,  2.8284,  2.8284]) 

在整个数组中找到标准误差,使用 n 自由度:

>>> stats.sem(a, axis=None, ddof=0)
1.2893796958227628 

scipy.stats.bayes_mvs

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

scipy.stats.bayes_mvs(data, alpha=0.9)

均值、方差和标准差的贝叶斯置信区间。

参数:

dataarray_like

输入数据,如果是多维的,则通过bayes_mvs将其展平为 1-D。需要至少 2 个数据点。

alpha浮点数,可选

返回置信区间包含真实参数的概率。

返回:

mean_cntr, var_cntr, std_cntr元组

这三个结果分别是均值、方差和标准差的元组形式:

(center, (lower, upper)) 

对于center,是给定数据的条件概率密度函数均值,对于(lower, upper),是以中位数为中心的置信区间,包含到概率alpha的估计。

另见

mvsdist

注意事项

每个均值、方差和标准差估计的元组表示为(center, (lower, upper)),其中 center 是给定数据的条件概率密度函数均值,(lower, upper)是以中位数为中心的置信区间,包含到概率alpha的估计。

转换数据为 1-D,假设所有数据具有相同的均值和方差。使用杰弗里先验法进行方差和标准差估计。

等效于tuple((x.mean(), x.interval(alpha)) for x in mvsdist(dat))

参考文献

T.E. Oliphant, “从数据中估计均值、方差和标准差的贝叶斯视角”,scholarsarchive.byu.edu/facpub/278,2006 年。

示例

首先是一个基本示例,用于展示输出:

>>> from scipy import stats
>>> data = [6, 9, 12, 7, 8, 8, 13]
>>> mean, var, std = stats.bayes_mvs(data)
>>> mean
Mean(statistic=9.0, minmax=(7.103650222612533, 10.896349777387467))
>>> var
Variance(statistic=10.0, minmax=(3.176724206..., 24.45910382...))
>>> std
Std_dev(statistic=2.9724954732045084,
 minmax=(1.7823367265645143, 4.945614605014631)) 

现在我们生成一些正态分布的随机数据,并使用 95%置信区间对均值和标准差的估计进行如下操作:

>>> n_samples = 100000
>>> data = stats.norm.rvs(size=n_samples)
>>> res_mean, res_var, res_std = stats.bayes_mvs(data, alpha=0.95) 
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.hist(data, bins=100, density=True, label='Histogram of data')
>>> ax.vlines(res_mean.statistic, 0, 0.5, colors='r', label='Estimated mean')
>>> ax.axvspan(res_mean.minmax[0],res_mean.minmax[1], facecolor='r',
...            alpha=0.2, label=r'Estimated mean (95% limits)')
>>> ax.vlines(res_std.statistic, 0, 0.5, colors='g', label='Estimated scale')
>>> ax.axvspan(res_std.minmax[0],res_std.minmax[1], facecolor='g', alpha=0.2,
...            label=r'Estimated scale (95% limits)') 
>>> ax.legend(fontsize=10)
>>> ax.set_xlim([-4, 4])
>>> ax.set_ylim([0, 0.5])
>>> plt.show() 

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

scipy.stats.mvsdist

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

scipy.stats.mvsdist(data)

数据的“冻结”分布,包括均值、方差和标准差。

参数:

dataarray_like

输入数组。使用 ravel 转换为 1-D。需要 2 个或更多数据点。

返回:

mdist“冻结”分布对象

表示数据均值的分布对象。

vdist“冻结”分布对象

表示数据方差的分布对象。

sdist“冻结”分布对象

表示数据标准差的分布对象。

参见

bayes_mvs

注意

bayes_mvs(data) 的返回值等同于 tuple((x.mean(), x.interval(0.90)) for x in mvsdist(data))

换句话说,在从此函数返回的三个分布对象上调用 <dist>.mean()<dist>.interval(0.90) 将返回与 bayes_mvs 返回的相同结果。

参考文献

T.E. Oliphant,“从数据中估计均值、方差和标准差的贝叶斯视角”,scholarsarchive.byu.edu/facpub/278,2006 年。

示例

>>> from scipy import stats
>>> data = [6, 9, 12, 7, 8, 8, 13]
>>> mean, var, std = stats.mvsdist(data) 

现在我们有了冻结的分布对象“mean”、“var”和“std”,我们可以进行检查:

>>> mean.mean()
9.0
>>> mean.interval(0.95)
(6.6120585482655692, 11.387941451734431)
>>> mean.std()
1.1952286093343936 

scipy.stats.entropy

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

scipy.stats.entropy(pk, qk=None, base=None, axis=0, *, nan_policy='propagate', keepdims=False)

计算给定分布的 Shannon 熵/相对熵。

如果仅提供了概率pk,则香农熵计算为H = -sum(pk * log(pk))

如果qk不为 None,则计算相对熵D = sum(pk * log(pk / qk))。这个量也被称为 Kullback-Leibler 散度。

如果pkqk的和不为 1,则此例程将对它们进行标准化。

参数:

pkarray_like

定义(离散)分布。对于pk的每个轴切片,元素i是事件i的(可能未标准化的)概率。

qkarray_like,可选

用于计算相对熵的序列。应与pk具有相同的格式。

basefloat,可选

要使用的对数基数,默认为e(自然对数)。

axisint 或 None,默认值:0

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

nan_policy

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

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

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

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

keepdimsbool,默认值:False

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

返回值:

S

计算得到的熵。

注意事项

通俗地讲,香农熵量化了离散随机变量可能结果的预期不确定性。例如,如果要对由一组符号序列组成的消息进行编码并通过无噪声信道传输,则香农熵H(pk)给出了每个符号所需的信息单位数的平均下界,如果符号的发生频率由离散分布pk控制[1]。基数的选择确定了单位的选择;例如,自然对数e用于 nats,2用于 bits,等等。

相对熵 D(pk|qk) 量化了如果编码针对概率分布 qk 而不是真实分布 pk 进行了优化,则每个符号所需的平均信息单位数的增加量。非正式地,相对熵量化了在真实分布实际为 pk 时,但人们认为其为 qk 时所经历的预期惊讶的过量。

相关量,交叉熵 CE(pk, qk),满足方程 CE(pk, qk) = H(pk) + D(pk|qk),也可以用公式 CE = -sum(pk * log(qk)) 计算。如果编码针对概率分布 qk 进行了优化,当真实分布为 pk 时,它给出每个符号所需的平均信息单位数。它不是直接由 entropy 计算的,但可以通过两次调用函数来计算(见示例)。

更多信息请参见 [2]

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

参考文献

[1]

Shannon, C.E. (1948),A Mathematical Theory of Communication. Bell System Technical Journal, 27: 379-423. doi.org/10.1002/j.1538-7305.1948.tb01338.x

[2]

Thomas M. Cover 和 Joy A. Thomas. 2006. Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing). Wiley-Interscience, USA.

示例

公平硬币的结果是最不确定的:

>>> import numpy as np
>>> from scipy.stats import entropy
>>> base = 2  # work in units of bits
>>> pk = np.array([1/2, 1/2])  # fair coin
>>> H = entropy(pk, base=base)
>>> H
1.0
>>> H == -np.sum(pk * np.log(pk)) / np.log(base)
True 

有偏硬币的结果不那么不确定:

>>> qk = np.array([9/10, 1/10])  # biased coin
>>> entropy(qk, base=base)
0.46899559358928117 

公平硬币和有偏硬币之间的相对熵计算如下:

>>> D = entropy(pk, qk, base=base)
>>> D
0.7369655941662062
>>> D == np.sum(pk * np.log(pk/qk)) / np.log(base)
True 

交叉熵可以计算为熵和相对熵的总和`:

>>> CE = entropy(pk, base=base) + entropy(pk, qk, base=base)
>>> CE
1.736965594166206
>>> CE == -np.sum(pk * np.log(qk)) / np.log(base)
True 

scipy.stats.differential_entropy

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

scipy.stats.differential_entropy(values, *, window_length=None, base=None, axis=0, method='auto', nan_policy='propagate', keepdims=False)

给定分布的样本,估计微分熵。

根据样本大小选择默认的方法,可使用method参数选择多种估计方法。

参数:

values序列

从连续分布中抽取样本。

window_length整数,可选

用于计算 Vasicek 估计的窗口长度。必须是 1 到样本大小的一半之间的整数。如果为None(默认值),则使用启发式值

[\left \lfloor \sqrt{n} + 0.5 \right \rfloor]

其中 (n) 是样本大小。这一启发式方法最初是在文献中提出的[2],现在在文献中很常见。

base浮点数,可选

使用的对数基数,默认为e(自然对数)。

axisint 或 None,默认值为:0

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

method,可选

从样本中估计微分熵的方法。默认为'auto'。更多信息请参阅注意事项。

nan_policy

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

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

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

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

keepdimsbool,默认值为:False

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

返回:

entropy浮点数

计算得到的微分熵。

注意事项

该函数在极限情况下将收敛到真实的微分熵

[n \to \infty, \quad m \to \infty, \quad \frac{m}{n} \to 0]

对于给定的样本大小,window_length的最佳选择取决于(未知的)分布。通常,分布的密度越平滑,window_length的最佳值就越大[1]

method参数有以下选项可供选择。

  • 'vasicek'使用[1]中提出的估计器。这是微分熵的最早和最有影响力的估计器之一。

  • 'van es'使用在[3]中提出的修正偏差估计器,不仅是一致的,而且在某些条件下渐近正态。

  • 'ebrahimi'使用在[4]中提出的估计器,在模拟中显示比 Vasicek 估计器具有更小的偏差和均方误差。

  • 'correa'使用在[5]中基于局部线性回归提出的估计器。在模拟研究中,其均方误差始终比 Vasiceck 估计器小,但计算成本更高。

  • 'auto'自动选择方法(默认)。目前,这为非常小的样本(<10)选择'van es',对于中等样本大小(11-1000)选择'ebrahimi',对于较大样本选择'vasicek',但此行为可能在未来版本中更改。

所有估计器均按照[6]中描述的方式实现。

从 SciPy 1.9 开始,不推荐新代码使用np.matrix输入,在进行计算之前将其转换为np.ndarray。在这种情况下,输出将是适当形状的标量或np.ndarray,而不是 2D 的np.matrix。类似地,尽管忽略了掩码数组的掩码元素,输出将是标量或np.ndarray,而不是带有mask=False的掩码数组。

参考文献

[1] (1,2)

Vasicek, O. (1976). 基于样本熵的正态性检验. 《皇家统计学会杂志:B 系列(方法学)》,38(1),54-59。

[2]

Crzcgorzewski, P., & Wirczorkowski, R. (1999). 基于熵的指数分布适合性检验. 《统计学通信-理论与方法》,28(5),1183-1202。

[3]

Van Es, B. (1992). 通过基于间隔的统计量类估计密度相关的函数. 《斯堪的纳维亚统计学杂志》,61-72。

[4]

Ebrahimi, N., Pflughoeft, K., & Soofi, E. S. (1994). 两种样本熵测量. 《统计与概率信函》,20(3),225-234。

[5]

Correa, J. C. (1995). 新的熵估计器. 《统计学通信-理论与方法》,24(10),2439-2449。

[6]

Noughabi, H. A. (2015). 使用数值方法进行熵估计. 《数据科学年鉴》,2(2),231-241。link.springer.com/article/10.1007/s40745-015-0045-9

示例

>>> import numpy as np
>>> from scipy.stats import differential_entropy, norm 

标准正态分布的熵:

>>> rng = np.random.default_rng()
>>> values = rng.standard_normal(100)
>>> differential_entropy(values)
1.3407817436640392 

与真实熵比较:

>>> float(norm.entropy())
1.4189385332046727 

对于在 5 到 1000 之间的多个样本大小,比较'vasicek''van es''ebrahimi'方法的准确性。具体比较(1000 次试验中)估计与分布真实差分熵之间的均方根误差。

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>>
>>>
>>> def rmse(res, expected):
...  '''Root mean squared error'''
...     return np.sqrt(np.mean((res - expected)**2))
>>>
>>>
>>> a, b = np.log10(5), np.log10(1000)
>>> ns = np.round(np.logspace(a, b, 10)).astype(int)
>>> reps = 1000  # number of repetitions for each sample size
>>> expected = stats.expon.entropy()
>>>
>>> method_errors = {'vasicek': [], 'van es': [], 'ebrahimi': []}
>>> for method in method_errors:
...     for n in ns:
...        rvs = stats.expon.rvs(size=(reps, n), random_state=rng)
...        res = stats.differential_entropy(rvs, method=method, axis=-1)
...        error = rmse(res, expected)
...        method_errors[method].append(error)
>>>
>>> for method, errors in method_errors.items():
...     plt.loglog(ns, errors, label=method)
>>>
>>> plt.legend()
>>> plt.xlabel('sample size')
>>> plt.ylabel('RMSE (1000 trials)')
>>> plt.title('Entropy Estimator Error (Exponential Distribution)') 

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

scipy.stats.median_abs_deviation

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

scipy.stats.median_abs_deviation(x, axis=0, center=<function median>, scale=1.0, nan_policy='propagate')

计算给定轴上数据的中位绝对偏差。

中位数绝对偏差(MAD,[1])计算从中位数到绝对偏差的中位数。这是一种与标准偏差类似但更鲁棒于异常值的离散度测量方法[2]

空数组的 MAD 是 np.nan

新版本 1.5.0 中新增。

参数:

x类似数组

可转换为数组的输入数组或对象。

axisint 或 None,可选

计算范围的轴。默认为 0。如果为 None,则在整个数组上计算 MAD。

中心可调用,可选

将返回中心值的函数。默认使用 np.median。任何用户定义的函数都需要具有 func(arr, axis) 的函数签名。

尺度标量或字符串,可选

尺度的数值将从最终结果中除去。默认为 1.0。还接受字符串“normal”,这将导致 scale 成为标准正态分位函数在 0.75 处的倒数,约为 0.67449。还允许类似数组的尺度,只要它正确广播到输出,使得 out / scale 是有效操作即可。输出维度取决于输入数组 xaxis 参数。

nan_policy,可选

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

  • ‘propagate’: 返回 nan。

  • ‘raise’: 抛出错误

  • ‘omit’: 在计算时忽略 nan 值。

返回:

mad标量或 ndarray

如果 axis=None,则返回一个标量。如果输入包含小于 np.float64 的整数或浮点数,则输出数据类型为 np.float64。否则,输出数据类型与输入相同。

另请参见

numpy.std, numpy.var, numpy.median, scipy.stats.iqr, scipy.stats.tmean

scipy.stats.tstd, scipy.stats.tvar

注释

center参数仅影响计算 MAD 时计算的中心值。也就是说,传入center=np.mean将计算围绕平均值的 MAD - 而不是计算平均绝对偏差。

输入数组可能包含inf,但如果center返回inf,则该数据对应的 MAD 将为nan

参考文献

[1]

“中位数绝对偏差”,en.wikipedia.org/wiki/Median_absolute_deviation

[2]

“尺度鲁棒性测量”,en.wikipedia.org/wiki/Robust_measures_of_scale

示例

当比较median_abs_deviationnp.std的行为时,后者在我们将数组的单个值更改为异常值时受影响,而 MAD 几乎没有变化:

>>> import numpy as np
>>> from scipy import stats
>>> x = stats.norm.rvs(size=100, scale=1, random_state=123456)
>>> x.std()
0.9973906394005013
>>> stats.median_abs_deviation(x)
0.82832610097857
>>> x[0] = 345.6
>>> x.std()
34.42304872314415
>>> stats.median_abs_deviation(x)
0.8323442311590675 

轴处理示例:

>>> x = np.array([[10, 7, 4], [3, 2, 1]])
>>> x
array([[10,  7,  4],
 [ 3,  2,  1]])
>>> stats.median_abs_deviation(x)
array([3.5, 2.5, 1.5])
>>> stats.median_abs_deviation(x, axis=None)
2.0 

标准化尺度示例:

>>> x = stats.norm.rvs(size=1000000, scale=2, random_state=123456)
>>> stats.median_abs_deviation(x)
1.3487398527041636
>>> stats.median_abs_deviation(x, scale='normal')
1.9996446978061115 

scipy.stats.cumfreq

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

scipy.stats.cumfreq(a, numbins=10, defaultreallimits=None, weights=None)

返回一个累积频率直方图,使用直方图函数。

累积直方图是一种映射,它计算了到指定箱子的所有箱子中的观测累积数。

参数:

aarray_like

输入数组。

numbinsint,可选

用于直方图的箱子数。默认为 10。

defaultreallimitstuple (lower, upper),可选

直方图的范围的下限和上限值。如果未指定值,则使用稍大于 a 值范围的范围。具体而言,(a.min() - s, a.max() + s),其中 s = (1/2)(a.max() - a.min()) / (numbins - 1)

weightsarray_like,可选

a 中每个值的权重。默认为 None,即每个值的权重为 1.0。

返回:

cumcountndarray

累积频率的分箱值。

lowerlimitfloat

较低的实际限制

binsizefloat

每个箱子的宽度。

extrapointsint

额外点。

示例

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> x = [1, 4, 2, 1, 3, 1]
>>> res = stats.cumfreq(x, numbins=4, defaultreallimits=(1.5, 5))
>>> res.cumcount
array([ 1.,  2.,  3.,  3.])
>>> res.extrapoints
3 

创建具有 1000 个随机值的正态分布

>>> samples = stats.norm.rvs(size=1000, random_state=rng) 

计算累积频率

>>> res = stats.cumfreq(samples, numbins=25) 

计算 x 的值的空间

>>> x = res.lowerlimit + np.linspace(0, res.binsize*res.cumcount.size,
...                                  res.cumcount.size) 

绘制直方图和累积直方图

>>> fig = plt.figure(figsize=(10, 4))
>>> ax1 = fig.add_subplot(1, 2, 1)
>>> ax2 = fig.add_subplot(1, 2, 2)
>>> ax1.hist(samples, bins=25)
>>> ax1.set_title('Histogram')
>>> ax2.bar(x, res.cumcount, width=res.binsize)
>>> ax2.set_title('Cumulative histogram')
>>> ax2.set_xlim([x.min(), x.max()]) 
>>> plt.show() 

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

scipy.stats.percentileofscore

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

scipy.stats.percentileofscore(a, score, kind='rank', nan_policy='propagate')

计算相对于一组分数的分数的百分位数。

例如,percentileofscore的 80%表示a中 80%的分数低于给定的分数。在存在间隙或并列值的情况下,确切的定义取决于可选关键字kind

参数:

aarray_like

与* score *进行比较的数组。

scorearray_like

用于计算百分位数的分数。

kind, optional

指定结果分数的解释。可用以下选项(默认为‘rank’):

  • ‘rank’:分数的平均百分位排名。在存在多个匹配项的情况下,平均所有匹配分数的百分位排名。
  • ‘weak’:这种类型对应于累积分布函数的定义。百分位数为 80%表示 80%的值小于或等于提供的分数。
  • ‘strict’:类似于“weak”,但仅计数严格小于给定分数的值。
  • ‘mean’:弱和严格分数的平均值,通常用于测试。参见 en.wikipedia.org/wiki/Percentile_rank

nan_policy, optional

指定如何处理a中的* nan *值。可用以下选项(默认为‘propagate’):

  • ‘propagate’:对于* score *中的每个值都返回 nan。
  • ‘raise’:抛出错误
  • ‘omit’:执行计算时忽略 nan 值

返回:

pcosfloat

分数在* a *中的百分位数(0-100)

另请参阅

numpy.percentile

scipy.stats.scoreatpercentile, scipy.stats.rankdata

示例

给定分数以下的给定值四分之三:

>>> import numpy as np
>>> from scipy import stats
>>> stats.percentileofscore([1, 2, 3, 4], 3)
75.0 

对于多个匹配项,请注意两个匹配项的分数分别为 0.6 和 0.8:

>>> stats.percentileofscore([1, 2, 3, 3, 4], 3)
70.0 

仅有 2/5 的值严格小于 3:

>>> stats.percentileofscore([1, 2, 3, 3, 4], 3, kind='strict')
40.0 

但是 4/5 的值小于或等于 3:

>>> stats.percentileofscore([1, 2, 3, 3, 4], 3, kind='weak')
80.0 

严格和弱得分之间的平均值是:

>>> stats.percentileofscore([1, 2, 3, 3, 4], 3, kind='mean')
60.0 

支持任意维度的分数数组:

>>> stats.percentileofscore([1, 2, 3, 3, 4], [2, 3])
array([40., 70.]) 

输入可以是无限的:

>>> stats.percentileofscore([-np.inf, 0, 1, np.inf], [1, 2, np.inf])
array([75., 75., 100.]) 

如果* a 为空,则生成的百分位数均为 nan *:

>>> stats.percentileofscore([], [1, 2])
array([nan, nan]) 

scipy.stats.scoreatpercentile

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

scipy.stats.scoreatpercentile(a, per, limit=(), interpolation_method='fraction', axis=None)

计算输入序列给定百分位数处的分数。

例如,per=50 处的分数是中位数。如果所需分位数位于两个数据点之间,我们根据 interpolation 的值进行插值。如果提供了 limit 参数,则应为两个值(下限,上限)的元组。

参数:

a类似数组

要提取分数的一维值数组。

per类似数组

要提取分数的百分位数。值应该在区间 [0,100] 内。

limit元组,可选

两个标量的元组,用于计算百分位数的下限和上限。a 的值如果在此(闭合)区间之外将被忽略。

interpolation_method,可选

指定在所需分位数位于两个数据点 ij 之间时要使用的插值方法。以下选项可用(默认为 ‘fraction’):

  • fraction: i + (j - i) * fraction,其中fraction是被ij包围的索引的分数部分
  • ‘lower’: i
  • ‘higher’: j

axis整数,可选

计算百分位数的轴。默认为 None。如果为 None,则在整个数组 a 上计算。

返回:

score浮点数或者 ndarray

百分位数处的分数。

另请参见

percentileofscorenumpy.percentile

注意

这个函数将来将会过时。对于 NumPy 1.9 及更高版本,建议使用 numpy.percentile,它提供了与 scoreatpercentile 相同的功能,并且速度显著更快。

示例

>>> import numpy as np
>>> from scipy import stats
>>> a = np.arange(100)
>>> stats.scoreatpercentile(a, 50)
49.5 

scipy.stats.relfreq

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

scipy.stats.relfreq(a, numbins=10, defaultreallimits=None, weights=None)

返回一个相对频率直方图,使用直方图函数。

相对频率直方图是每个区间内观察值数量相对于总观察值的映射。

参数:

a数组型

输入数组。

numbins整数,可选

直方图使用的箱子数量。默认为 10。

defaultreallimits元组(下限,上限),可选

直方图的范围的下限和上限值。如果未给定值,则使用稍大于 a 值范围的范围。具体来说是 (a.min() - s, a.max() + s),其中 s = (1/2)(a.max() - a.min()) / (numbins - 1)

weights数组型,可选

a中每个值的权重。默认为 None,每个值的权重为 1.0。

返回:

frequencyn 维数组

相对频率的分箱值。

lowerlimit浮点数

较低的实际限制。

binsize浮点数

每个箱子的宽度。

extrapoints整数

额外的点。

示例

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> a = np.array([2, 4, 1, 2, 3, 2])
>>> res = stats.relfreq(a, numbins=4)
>>> res.frequency
array([ 0.16666667, 0.5       , 0.16666667,  0.16666667])
>>> np.sum(res.frequency)  # relative frequencies should add up to 1
1.0 

创建具有 1000 个随机值的正态分布

>>> samples = stats.norm.rvs(size=1000, random_state=rng) 

计算相对频率

>>> res = stats.relfreq(samples, numbins=25) 

计算 x 的值空间

>>> x = res.lowerlimit + np.linspace(0, res.binsize*res.frequency.size,
...                                  res.frequency.size) 

绘制相对频率直方图

>>> fig = plt.figure(figsize=(5, 4))
>>> ax = fig.add_subplot(1, 1, 1)
>>> ax.bar(x, res.frequency, width=res.binsize)
>>> ax.set_title('Relative frequency histogram')
>>> ax.set_xlim([x.min(), x.max()]) 
>>> plt.show() 

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

scipy.stats.binned_statistic

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

scipy.stats.binned_statistic(x, values, statistic='mean', bins=10, range=None)

计算一个或多个数据集的分 bin 统计量。

这是直方图函数的一般化。直方图将空间划分为 bins,并返回每个 bin 中点的数量。此函数允许计算每个 bin 中的值(或值集合)的和、均值、中位数或其他统计量。

参数:

x(N,) 类似数组

要分 bin 的一系列值。

values(N,) 类似数组或者(N,) 类似数组的列表

统计量将被计算的数据。这必须与x具有相同的形状,或者是一组序列 - 每个序列与x具有相同的形状。如果values是一组序列,则将独立地计算每个统计量。

statistic字符串或可调用对象,可选

要计算的统计量(默认为‘mean’)。以下统计量可用:

  • ‘mean’:计算每个 bin 内点的平均值。空的 bins 将用 NaN 表示。
  • ‘std’:计算每个 bin 内的标准差。这是使用 ddof=0 隐式计算的。
  • ‘median’:计算每个 bin 内点的值的中位数。空的 bins 将用 NaN 表示。
  • ‘count’:计算每个 bin 内点的数量。这等同于一个非加权直方图。values数组不被引用。
  • ‘sum’:计算每个 bin 内点的值的总和。这等同于一个加权直方图。
  • ‘min’:计算每个 bin 内点的最小值。空的 bins 将用 NaN 表示。
  • ‘max’:计算每个 bin 内点的值的最大值。空的 bins 将用 NaN 表示。
  • function:一个用户定义的函数,接受一个值的 1D 数组,并输出一个单一的数值统计量。此函数将在每个 bin 中的值上调用。空的 bins 将由 function([])表示,如果这导致错误,则返回 NaN。

bins整数或标量序列,可选

如果bins是整数,则定义给定范围内的等宽 bin 数(默认为 10)。如果bins是序列,则定义 bin 边缘,包括右边的边缘,允许非均匀的 bin 宽度。小于最低 bin 边缘的x值被分配给 bin 号 0,超出最高 bin 的值被分配给bins[-1]。如果指定了 bin 边缘,则 bins 的数量将为(nx = len(bins)-1)。

range (float, float) 或 [(float, float)],可选

bins 的下限和上限。如果未提供,则范围为(x.min(), x.max())。超出范围的值将被忽略。

返回值:

statistic数组

每个 bin 中所选统计量的值。

bin_edges浮点数 dtype 的数组

返回 bin 边界 (length(statistic)+1)

binnumber:整数型的 1-D ndarray

每个值 x 属于的箱子(对应于 bin_edges)的索引。与 values 长度相同。箱号为 i 表示对应的值位于 (bin_edges[i-1], bin_edges[i]) 之间。

另请参阅

numpy.digitize, numpy.histogram, binned_statistic_2d, binned_statistic_dd

注意:

除了最后一个(最右边的)箱子是半开放的。换句话说,如果 bins[1, 2, 3, 4],那么第一个箱子是 [1, 2)(包括 1,但不包括 2),第二个是 [2, 3)。然而,最后一个箱子是 [3, 4],其中包括 4。

版本 0.11.0 中的新特性。

示例:

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

首先是一些基本示例:

在给定样本范围内创建两个均匀间隔的箱子,并计算每个箱子中对应的值的总和:

>>> values = [1.0, 1.0, 2.0, 1.5, 3.0]
>>> stats.binned_statistic([1, 1, 2, 5, 7], values, 'sum', bins=2)
BinnedStatisticResult(statistic=array([4\. , 4.5]),
 bin_edges=array([1., 4., 7.]), binnumber=array([1, 1, 1, 2, 2])) 

也可以传递多个值数组。统计量是在每个集合上独立计算的:

>>> values = [[1.0, 1.0, 2.0, 1.5, 3.0], [2.0, 2.0, 4.0, 3.0, 6.0]]
>>> stats.binned_statistic([1, 1, 2, 5, 7], values, 'sum', bins=2)
BinnedStatisticResult(statistic=array([[4\. , 4.5],
 [8\. , 9\. ]]), bin_edges=array([1., 4., 7.]),
 binnumber=array([1, 1, 1, 2, 2])) 
>>> stats.binned_statistic([1, 2, 1, 2, 4], np.arange(5), statistic='mean',
...                        bins=3)
BinnedStatisticResult(statistic=array([1., 2., 4.]),
 bin_edges=array([1., 2., 3., 4.]),
 binnumber=array([1, 2, 1, 2, 3])) 

作为第二个例子,我们现在生成一些作为风速函数的帆船速度的随机数据,然后确定我们的船在特定风速下的速度有多快:

>>> rng = np.random.default_rng()
>>> windspeed = 8 * rng.random(500)
>>> boatspeed = .3 * windspeed**.5 + .2 * rng.random(500)
>>> bin_means, bin_edges, binnumber = stats.binned_statistic(windspeed,
...                 boatspeed, statistic='median', bins=[1,2,3,4,5,6,7])
>>> plt.figure()
>>> plt.plot(windspeed, boatspeed, 'b.', label='raw data')
>>> plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='g', lw=5,
...            label='binned statistic of data')
>>> plt.legend() 

现在我们可以使用 binnumber 来选择所有风速低于 1 的数据点:

>>> low_boatspeed = boatspeed[binnumber == 0] 

最后一个例子中,我们将使用 bin_edgesbinnumber 来绘制一个分布的图,该图显示每个 bin 中的平均值及其周围的分布,叠加在常规直方图和概率分布函数之上:

>>> x = np.linspace(0, 5, num=500)
>>> x_pdf = stats.maxwell.pdf(x)
>>> samples = stats.maxwell.rvs(size=10000) 
>>> bin_means, bin_edges, binnumber = stats.binned_statistic(x, x_pdf,
...         statistic='mean', bins=25)
>>> bin_width = (bin_edges[1] - bin_edges[0])
>>> bin_centers = bin_edges[1:] - bin_width/2 
>>> plt.figure()
>>> plt.hist(samples, bins=50, density=True, histtype='stepfilled',
...          alpha=0.2, label='histogram of data')
>>> plt.plot(x, x_pdf, 'r-', label='analytical pdf')
>>> plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='g', lw=2,
...            label='binned statistic of data')
>>> plt.plot((binnumber - 0.5) * bin_width, x_pdf, 'g.', alpha=0.5)
>>> plt.legend(fontsize=10)
>>> plt.show() 

../../_images/scipy-stats-binned_statistic-1_00.png../../_images/scipy-stats-binned_statistic-1_01.png

scipy.stats.binned_statistic_2d

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

scipy.stats.binned_statistic_2d(x, y, values, statistic='mean', bins=10, range=None, expand_binnumbers=False)

计算一个或多个数据集的二维分箱统计量。

这是直方图 2d 函数的泛化。直方图将空间划分为箱子,并返回每个箱子中点的数量。此函数允许计算每个箱内值(或值集)的总和、均值、中位数或其他统计量。

参数:

x(N,) 数组样式

要沿第一维度分箱的一系列值。

y(N,) 数组样式

要沿第二维度分箱的一系列值。

values(N,) 数组样式或(N,)数组样式的列表

将计算统计量的数据。它必须与x的形状相同,或者是一个序列列表 - 每个序列具有与x相同的形状。如果values是这样的列表,则将分别在每个列表上计算统计量。

statistic字符串或可调用对象,可选

要计算的统计量(默认为‘mean’)。可用以下统计量:

  • ‘mean’:计算每个箱内点的平均值。空箱将用 NaN 表示。
  • ‘std’:计算每个箱内的标准偏差。这隐式地以 ddof=0 计算。
  • ‘median’:计算每个箱内点的中位数。空箱将用 NaN 表示。
  • ‘count’:计算每个箱内的点数。这与未加权的直方图相同。values数组未被引用。
  • ‘sum’:计算每个箱内点的总和。这与加权直方图相同。
  • ‘min’:计算每个箱内点的最小值。空箱将用 NaN 表示。
  • ‘max’:计算每个箱内点的最大值。空箱将用 NaN 表示。
  • 函数:一个用户定义的函数,接受值的 1D 数组,并输出单个数值统计量。该函数将在每个箱子的值上调用。空箱将由 function([])表示,或者如果返回错误,则为 NaN。

bins整数或[int, int]或类数组或[array, array],可选

箱子的规范:

  • 两个维度的箱子数量(nx = ny = bins),
  • 每个维度中的箱子数量(nx, ny = bins),
  • 两个维度的箱子边缘(x_edge = y_edge = bins),
  • 每个维度中的箱子边缘(x_edge, y_edge = bins)。

如果指定了箱子边缘,则箱子数量将为(nx = len(x_edge)-1, ny = len(y_edge)-1)。

range(2,2) 数组样式,可选

每个维度的最左边和最右边的箱子边缘(如果在bins参数中未显式指定):[[xmin, xmax], [ymin, ymax]]。该范围外的所有值将被视为异常值,并不计入直方图。

expand_binnumbers布尔值,可选

'False'(默认值):返回的 binnumber 是形状为 (N,) 的线性化 bin 索引数组。'True':返回的 binnumber '展开' 成形状为 (2,N) 的 ndarray,其中每行给出相应维度中的 bin 数。参见 binnumber 返回值和 Examples 部分。

新版本 0.17.0 中新增。

返回:

statistic(nx, ny) ndarray

每个二维 bin 中所选统计量的值。

x_edge(nx + 1) ndarray

第一维度上的 bin 边缘。

y_edge(ny + 1) ndarray

第二维度上的 bin 边缘。

binnumber(N,) 整数数组或 (2,N) 整数 ndarray

此参数为每个 sample 元素分配一个整数,表示此观察值所在的 bin。表示依赖于 expand_binnumbers 参数。详细信息请参见 Notes

另请参阅

numpy.digitizenumpy.histogram2dbinned_statisticbinned_statistic_dd

Notes

Binedges:除了最后一个(最右侧)bin 是半开放的。换句话说,如果 bins[1, 2, 3, 4],则第一个 bin 是 [1, 2)(包括 1,但不包括 2),第二个是 [2, 3)。然而,最后一个 bin 是 [3, 4],包含 4。

binnumber:此返回参数为每个 sample 元素分配一个整数,表示其所属的 bin。表示依赖于 expand_binnumbers 参数。如果为 'False'(默认值):返回的 binnumber 是一个形状为 (N,) 的数组,其中包含将每个 sample 元素映射到其对应 bin 的线性化索引(使用行优先顺序)。请注意,返回的线性化 bin 索引用于具有外部 bin 边界上额外 bin 的数组,以捕获定义的 bin 边界之外的值。如果为 'True':返回的 binnumber 是一个形状为 (2,N) 的 ndarray,其中每行分别指示每个维度的 bin 放置。在每个维度中,binnumber 为 i 表示相应的值位于(D_edge[i-1],D_edge[i])之间,其中 'D' 可以是 'x' 或 'y'。

新版本 0.11.0 中新增。

示例

>>> from scipy import stats 

使用显式 bin 边缘计算计数:

>>> x = [0.1, 0.1, 0.1, 0.6]
>>> y = [2.1, 2.6, 2.1, 2.1]
>>> binx = [0.0, 0.5, 1.0]
>>> biny = [2.0, 2.5, 3.0]
>>> ret = stats.binned_statistic_2d(x, y, None, 'count', bins=[binx, biny])
>>> ret.statistic
array([[2., 1.],
 [1., 0.]]) 

每个样本所在的 bin 由返回的 binnumber 参数给出。默认情况下,这些是线性化的 bin 索引:

>>> ret.binnumber
array([5, 6, 5, 9]) 

也可以使用 expand_binnumbers 参数将 bin 索引扩展为每个维度的单独条目:

>>> ret = stats.binned_statistic_2d(x, y, None, 'count', bins=[binx, biny],
...                                 expand_binnumbers=True)
>>> ret.binnumber
array([[1, 1, 1, 2],
 [1, 2, 1, 1]]) 

表明前三个元素属于 xbin 1,第四个元素属于 xbin 2;y 以此类推。

scipy.stats.binned_statistic_dd

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

scipy.stats.binned_statistic_dd(sample, values, statistic='mean', bins=10, range=None, expand_binnumbers=False, binned_statistic_result=None)

计算一组数据的多维度分箱统计量。

这是 histogramdd 函数的泛化。直方图将空间划分为箱,并返回每个箱内点的计数。此函数允许计算每个箱内值的总和、平均值、中位数或其他统计量。

参数:

sample 数组样本

将作为 N 个长度为 D 的数组序列或者是(N,D)数组传递的数据直方图。

values (N,) 数组样本或 (N,) 数组样本列表

将计算统计量的数据。它必须与sample具有相同的形状,或者是一个序列的列表 - 每个序列都具有与sample相同的形状。如果values是这样的列表,则将独立地对每个计算统计量。

statistic 字符串或可调用对象,可选

要计算的统计量(默认为'mean')。提供以下统计量选项:

  • ‘mean’:计算每个箱内点的值的平均值。空箱将由 NaN 表示。
  • ‘median’:计算每个箱内点的值的中位数。空箱将由 NaN 表示。
  • ‘count’:计算每个箱内点的计数。这等同于一个非加权直方图。values 数组不被引用。
  • ‘sum’:计算每个箱内点的值的总和。这等同于加权直方图。
  • ‘std’:计算每个箱内的标准差。这是以 ddof=0 隐式计算的。如果给定箱内值的数量为 0 或 1,则计算的标准差值将为该箱的 0。
  • ‘min’:计算每个箱内点的值的最小值。空箱将由 NaN 表示。
  • ‘max’:计算每个箱内点的值的最大值。空箱将由 NaN 表示。
  • function:一个用户定义的函数,接受一个 1D 值数组,并输出单个数值统计量。将对每个箱中的值调用此函数。空箱将由 function([]) 表示,如果这导致错误,则为 NaN。

bins 序列或正整数,可选

箱规范必须是以下形式之一:

  • 一系列描述沿每个维度的箱边界的数组。
  • 每个维度的箱数(nx、ny,... = bins)。
  • 所有维度的箱数(nx = ny = ... = bins)。

range 序列,可选

一系列下限和上限箱边界,如果在bins中未显式给出边界,则默认为每个维度的最小值和最大值。

expand_binnumbers 布尔值,可选

‘False’(默认值):返回的 binnumber 是一个形状为 (N,) 的数组,其中包含线性化的箱索引。‘True’:返回的 binnumber 被‘展开’为一个形状为 (D,N) 的 ndarray,其中每行给出相应维度的箱编号。请参见 binned_statistic_2d 中的 binnumber 返回值和 Examples 部分。

binned_statistic_result:binnedStatisticddResult

前一次调用函数的结果,以便在新值和/或不同统计信息上重用箱边界和箱号。若要重用箱号,expand_binnumbers 必须设置为 False(默认值)。

从版本 0.17.0 开始新增。

返回:

statistic:ndarray,形状为 (nx1, nx2, nx3,…)

每个二维箱中所选统计量的值。

bin_edges 的列表,其中包含 ndarray

描述每个维度的 (nxi + 1) 箱边缘的 D 个数组的列表。

binnumber(N,) 整数数组或 (D,N) 整数数组

此分配给 sample 中每个元素一个整数,表示此观测值所在的箱。表示依赖于 expand_binnumbers 参数。有关详细信息,请参见 Notes

另请参阅

numpy.digitize, numpy.histogramdd, binned_statistic, binned_statistic_2d

注释

Binedges:除最后一个(最右边的)箱子外,在每个维度中都是半开放的。换句话说,如果 bins[1, 2, 3, 4],那么第一个箱是 [1, 2)(包括 1,但不包括 2),第二个箱是 [2, 3)。然而,最后一个箱是 [3, 4],其中包含 4。

binnumber:此返回的参数将为 sample 中的每个元素分配一个整数,表示其所属的箱。表示依赖于 expand_binnumbers 参数。若‘False’(默认值):返回的 binnumber 是一个形状为 (N,) 的数组,其中每个元素的索引映射到其相应的箱(使用行优先顺序)。若‘True’:返回的 binnumber 是一个形状为 (D,N) 的 ndarray,其中每行分别指示每个维度的箱位置。在每个维度中,一个 binnumberi 意味着对应的值在(bin_edges[D][i-1], bin_edges[D][i])之间。

从版本 0.11.0 开始新增。

示例

>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D 

以一个包含 600 个 (x, y) 坐标的数组为例。binned_statistic_dd 可以处理更高维度 D 的数组。但需要维度 D+1 的绘图。

>>> mu = np.array([0., 1.])
>>> sigma = np.array([[1., -0.5],[-0.5, 1.5]])
>>> multinormal = stats.multivariate_normal(mu, sigma)
>>> data = multinormal.rvs(size=600, random_state=235412)
>>> data.shape
(600, 2) 

创建箱子并计算每个箱子中有多少数组:

>>> N = 60
>>> x = np.linspace(-3, 3, N)
>>> y = np.linspace(-3, 4, N)
>>> ret = stats.binned_statistic_dd(data, np.arange(600), bins=[x, y],
...                                 statistic='count')
>>> bincounts = ret.statistic 

设置柱的体积和位置:

>>> dx = x[1] - x[0]
>>> dy = y[1] - y[0]
>>> x, y = np.meshgrid(x[:-1]+dx/2, y[:-1]+dy/2)
>>> z = 0 
>>> bincounts = bincounts.ravel()
>>> x = x.ravel()
>>> y = y.ravel() 
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111, projection='3d')
>>> with np.errstate(divide='ignore'):   # silence random axes3d warning
...     ax.bar3d(x, y, z, dx, dy, bincounts) 

重复使用具有新值的箱编号和箱边缘:

>>> ret2 = stats.binned_statistic_dd(data, -np.arange(600),
...                                  binned_statistic_result=ret,
...                                  statistic='mean') 

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

scipy.stats.ttest_1samp

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

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

计算一个组分数的均值的 T 检验。

这是一个关于期望值(样本的平均值)的空假设的检验,即样本 a 的期望值等于给定的总体均值 popmean

参数:

a类似数组

样本观察值。

popmeanfloat 或 类似数组

空假设中的期望值。如果是类似数组,则其沿 axis 的长度必须等于 1,否则必须可以广播至 a

axis整数或 None,默认为:0

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

nan_policy

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

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

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

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

alternative,可选

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

  • ‘two-sided’: 样本的基础分布的均值与给定的总体均值不同(popmean

  • ‘less’: 样本的基础分布的均值小于给定的总体均值(popmean

  • ‘greater’: 样本的基础分布的均值大于给定的总体均值(popmean

keepdims布尔值,默认为:False

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

返回:

resultTtestResult

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

statisticfloat 或 数组

t 统计量。

pvaluefloat 或 数组

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

dffloat 或 数组

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

在版本 1.10.0 中新增。

对象还具有以下方法:

confidence_interval(confidence_level=0.95)

计算给定置信水平下围绕总体均值的置信区间。置信区间以具有 lowhigh 字段的 namedtuple 返回。

在版本 1.10.0 中新增。

注意事项

统计量计算公式为(np.mean(a) - popmean)/se,其中se表示标准误差。因此,当样本均值大于总体均值时,统计量为正;当样本均值小于总体均值时,统计量为负。

从 SciPy 1.9 开始,np.matrix输入(不建议新代码使用)在执行计算之前会转换为np.ndarray。在这种情况下,输出将是一个标量或适当形状的np.ndarray,而不是 2D 的np.matrix。同样,虽然被屏蔽的数组的屏蔽元素会被忽略,但输出将是一个标量或np.ndarray,而不是带有mask=False的屏蔽数组。

例子

假设我们希望测试总体均值等于 0.5 的空假设。我们选择 99%的置信水平;也就是说,如果 p 值小于 0.01,我们将拒绝空假设,支持备选假设。

在进行来自标准均匀分布的随机变量测试时,该分布均值为 0.5,我们预期数据大多数时间与空假设一致。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> rvs = stats.uniform.rvs(size=50, random_state=rng)
>>> stats.ttest_1samp(rvs, popmean=0.5)
TtestResult(statistic=2.456308468440, pvalue=0.017628209047638, df=49) 

如预期的那样,0.017 的 p 值不低于我们的 0.01 阈值,因此我们不能拒绝空假设。

在测试来自标准正态分布的数据时,其均值为 0,我们预期将拒绝空假设。

>>> rvs = stats.norm.rvs(size=50, random_state=rng)
>>> stats.ttest_1samp(rvs, popmean=0.5)
TtestResult(statistic=-7.433605518875, pvalue=1.416760157221e-09, df=49) 

确实,p 值低于我们的 0.01 阈值,因此我们拒绝空假设,支持默认的“双侧”替代假设:总体均值不等于0.5。

然而,假设我们针对单侧替代检验空假设,即总体均值大于0.5。由于标准正态分布的均值小于 0.5,我们不会期望拒绝空假设。

>>> stats.ttest_1samp(rvs, popmean=0.5, alternative='greater')
TtestResult(statistic=-7.433605518875, pvalue=0.99999999929, df=49) 

毫不奇怪,由于 p 值大于我们的阈值,我们不会拒绝空假设。

注意,在使用 99%置信水平时,真空假设将被拒绝约 1%的时间。

>>> rvs = stats.uniform.rvs(size=(100, 50), random_state=rng)
>>> res = stats.ttest_1samp(rvs, popmean=0.5, axis=1)
>>> np.sum(res.pvalue < 0.01)
1 

实际上,即使以上所有 100 个样本均来自标准均匀分布,其总体均值确实为 0.5,我们也会错误地拒绝一个样本的空假设。

ttest_1samp还可以计算围绕总体均值的置信区间。

>>> rvs = stats.norm.rvs(size=50, random_state=rng)
>>> res = stats.ttest_1samp(rvs, popmean=0)
>>> ci = res.confidence_interval(confidence_level=0.95)
>>> ci
ConfidenceInterval(low=-0.3193887540880017, high=0.2898583388980972) 

95%置信区间的边界是参数popmean的最小和最大值,使得测试的 p 值为 0.05。

>>> res = stats.ttest_1samp(rvs, popmean=ci.low)
>>> np.testing.assert_allclose(res.pvalue, 0.05)
>>> res = stats.ttest_1samp(rvs, popmean=ci.high)
>>> np.testing.assert_allclose(res.pvalue, 0.05) 

在关于从样本抽取的总体的某些假设下,95%置信水平的置信区间预计在 95%的样本复制中包含真实总体均值。

>>> rvs = stats.norm.rvs(size=(50, 1000), loc=1, random_state=rng)
>>> res = stats.ttest_1samp(rvs, popmean=0)
>>> ci = res.confidence_interval()
>>> contains_pop_mean = (ci.low < 1) & (ci.high > 1)
>>> contains_pop_mean.sum()
953 

scipy.stats.binomtest

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

scipy.stats.binomtest(k, n, p=0.5, alternative='two-sided')

执行成功概率为 p 的检验。

二项式检验[1]是对伯努利实验中成功概率为p的原假设的检验。

可在许多统计学文本中找到检验的详细信息,例如[2]的 24.5 节。

参数:

k整数

成功次数。

n整数

试验次数。

p浮点数,可选

成功的假设概率,即预期的成功比例。该值必须在区间0 <= p <= 1内。默认值为p = 0.5

alternative,可选

指示备择假设。默认值为'两侧'。

返回:

resultBinomTestResult 实例

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

k 整数

成功次数(从binomtest中复制)

n 整数

试验次数(从binomtest中复制)

alternativestr

指示输入到binomtest中的备择假设。它将是'two-sided''greater''less'中的一个。

统计浮点数

成功比例的估计。

p 值浮点数

假设检验的 p 值。

该对象具有以下方法:

proportion_ci(confidence_level=0.95, method=’exact’):

计算statistic的置信区间。

注意事项

新功能版本 1.7.0。

参考文献

[1]

二项式检验,en.wikipedia.org/wiki/Binomial_test

[2]

Jerrold H. Zar,《生物统计分析》(第五版),Prentice Hall,Upper Saddle River,New Jersey USA(2010)

示例

>>> from scipy.stats import binomtest 

汽车制造商声称他们的汽车不安全的比例不超过 10%。检查了 15 辆汽车的安全性,发现 3 辆不安全。检验制造商的声明:

>>> result = binomtest(3, n=15, p=0.1, alternative='greater')
>>> result.pvalue
0.18406106910639114 

在 5%显著水平下,无法拒绝原假设,因为返回的 p 值大于 5%的临界值。

检验统计量等于估计比例,即简单地3/15

>>> result.statistic
0.2 

我们可以使用结果的proportion_ci()方法计算估计的置信区间:

>>> result.proportion_ci(confidence_level=0.95)
ConfidenceInterval(low=0.05684686759024681, high=1.0) 

scipy.stats.quantile_test

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

scipy.stats.quantile_test(x, *, q=0, p=0.5, alternative='two-sided')

执行分位数测试并计算分位数的置信区间。

此函数测试零假设,即 q 是样本 x 底层人口分布的分位数值。例如,默认参数下,它测试 x 底层分布的中位数是否为零。函数返回一个对象,包括测试统计量、p 值以及计算分位数置信区间的方法。

参数:

xarray_like

一维样本。

qfloat,默认值:0

假设的分位数值。

pfloat,默认值:0.5

与分位数相关联的概率;即人口比例小于 qp。必须严格在 0 和 1 之间。

alternative,可选

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

  • ‘two-sided’: 与概率 p 相关的分位数不是 q

  • ‘less’: 与概率 p 相关的分位数小于 q

  • ‘greater’: 与概率 p 相关的分位数大于 q

返回值:

resultQuantileTestResult

具有以下属性的对象:

statisticfloat

两种可能用于分位数测试的检验统计量之一。第一个检验统计量 T1x 中小于或等于假设分位数 q 的样本比例。第二个检验统计量 T2x 中严格小于假设分位数 q 的样本比例。

alternative = 'greater' 时,使用 T1 计算 p 值,statistic 设置为 T1

alternative = 'less' 时,使用 T2 计算 p 值,statistic 设置为 T2

alternative = 'two-sided' 时,考虑 T1T2,并使用导致最小 p 值的那个。

statistic_typeint

根据使用 T1T2 计算 p 值而定,为 12

pvaluefloat

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

该对象还具有以下方法:

confidence_interval(confidence_level=0.95)

计算与人口分位数相关联的概率 p 的置信区间。置信区间在 namedtuple 中返回,字段为 lowhigh。当观测数不足以计算所需置信度的置信区间时,值为 nan

注:

此测试及其计算置信区间的方法是非参数的。仅当观测值独立同分布时才有效。

测试遵循 Conover [1]。考虑了两个检验统计量。

T1:小于或等于 q 的观测值数量。

T1 = (x <= q).sum()

T2:严格小于 q 的观测值数量。

T2 = (x < q).sum()

使用两个检验统计量是必要的,以处理 x 是从离散或混合分布生成的可能性。

检验的零假设是:

H0:第 (p^{\mathrm{th}}) 个总体分位数是 q

并且每个检验统计量的零分布是 (\mathrm{binom}\left(n, p\right))。当 alternative='less' 时,备择假设是:

H1:第 (p^{\mathrm{th}}) 个总体分位数小于 q

而 p 值是二项随机变量小于或等于观测值 T1 的概率。

[Y \sim \mathrm{binom}\left(n, p\right)]

大于或等于观测值 T2

alternative='greater' 时,备择假设是:

H1:第 (p^{\mathrm{th}}) 个总体分位数大于 q

而 p 值是二项随机变量 Y 小于或等于观测值 T1 的概率。

alternative='two-sided' 时,备择假设是:

H1:q 不是第 (p^{\mathrm{th}}) 个总体分位数。

而 p 值是 'less''greater' 情况下 p 值的较小者的两倍。对于同一数据,这两个 p 值都可能超过 0.5,因此该值被限制在区间 ([0, 1])。

置信区间的方法归因于 Thompson [2],并且后来被证明适用于任何一组独立同分布样本 [3]。计算基于这样的观察:分位数 (q) 大于任何观测值 (x_m (1\leq m \leq N)) 的概率可以计算为

[\mathbb{P}(x_m \leq q) = 1 - \sum_{k=0}^{m-1} \binom{N}{k} qk(1-q)]

默认情况下,置信区间是针对 95% 的置信水平计算的。对于 95% 置信区间的常见解释是,如果从同一总体中重复抽取独立同分布样本并形成置信区间,这些置信区间将在大约 95% 的试验中包含指定分位数的真实值。

QuantileNPCI R 包中有类似的功能 [4]。其基础相同,但通过在样本值之间进行插值来计算置信区间边界,而本函数仅使用样本值作为边界。因此,quantile_test.confidence_interval 返回更保守的区间(即更大)。

计算分位数置信区间的相同方法包含在 confintr 包中[5]

双侧置信区间不能保证是最优的;即可能存在一个更紧的区间,其概率大于置信水平包含感兴趣的分位数。在没有关于样本的进一步假设(例如,底层分布的性质)的情况下,单侧区间是最优紧的。

参考文献

[1]

    1. Conover,《实用非参数统计学》,第 3 版,1999 年。

[2]

W. R. Thompson,《关于中位数和其他期望分布的置信区间》,《数理统计学年刊》,卷 7,第 3 期,pp. 122-128,1936 年,访问日期:2019 年 9 月 18 日。[在线]. 可用:www.jstor.org/stable/2957563.

[3]

H. A. David 和 H. N. Nagaraja,《非参数推断中的序统计量》,《序统计量》,John Wiley & Sons, Ltd,2005 年,pp. 159-170. 可用:onlinelibrary.wiley.com/doi/10.1002/0471722162.ch7.

[4]

N. Hutson, A. Hutson, L. Yan,《QuantileNPCI: 非参数置信区间的分位数》,R 包,cran.r-project.org/package=QuantileNPCI

[5]

M. Mayer,《confintr: 置信区间》,R 包,cran.r-project.org/package=confintr

例子

假设我们希望测试一个总体中位数等于 0.5 的零假设。我们选择 99%的置信水平;也就是说,如果 p 值小于 0.01,我们将拒绝零假设,接受备择假设。

当测试来自标准均匀分布的随机变量时,其中位数为 0.5,我们预期数据大部分时间与零假设一致。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> rvs = stats.uniform.rvs(size=100, random_state=rng)
>>> stats.quantile_test(rvs, q=0.5, p=0.5)
QuantileTestResult(statistic=45, statistic_type=1, pvalue=0.36820161732669576) 

顾名思义,p 值未低于我们的 0.01 阈值,因此我们无法拒绝零假设。

当测试来自标准正态分布的数据时,其中位数为 0,我们预期会拒绝零假设。

>>> rvs = stats.norm.rvs(size=100, random_state=rng)
>>> stats.quantile_test(rvs, q=0.5, p=0.5)
QuantileTestResult(statistic=67, statistic_type=2, pvalue=0.0008737198369123724) 

实际上,p 值低于我们的 0.01 阈值,因此我们拒绝零假设,接受默认的“双边”备择假设:总体中位数等于 0.5。

然而,假设我们要测试一个单侧备择假设,即总体中位数大于0.5。由于标准正态分布的中位数小于 0.5,我们不希望拒绝零假设。

>>> stats.quantile_test(rvs, q=0.5, p=0.5, alternative='greater')
QuantileTestResult(statistic=67, statistic_type=1, pvalue=0.9997956114162866) 

不出所料,由于 p 值高于我们的阈值,我们不会拒绝零假设,而是接受所选择的备择假设。

分位数测试不仅可以用于中位数,还可以用于任何分位数。例如,我们可以测试样本基础分布的第三四分位数是否大于 0.6。

>>> rvs = stats.uniform.rvs(size=100, random_state=rng)
>>> stats.quantile_test(rvs, q=0.6, p=0.75, alternative='greater')
QuantileTestResult(statistic=64, statistic_type=1, pvalue=0.00940696592998271) 

p 值低于阈值。我们拒绝零假设,接受备择假设:样本基础分布的第三四分位数大于 0.6。

quantile_test 还可以计算任何分位数的置信区间。

>>> rvs = stats.norm.rvs(size=100, random_state=rng)
>>> res = stats.quantile_test(rvs, q=0.6, p=0.75)
>>> ci = res.confidence_interval(confidence_level=0.95)
>>> ci
ConfidenceInterval(low=0.284491604437432, high=0.8912531024914844) 

在测试单侧备择假设时,置信区间包含所有观察结果,使得如果作为q,则测试的 p 值大于 0.05,因此不会拒绝原假设。例如:

>>> rvs.sort()
>>> q, p, alpha = 0.6, 0.75, 0.95
>>> res = stats.quantile_test(rvs, q=q, p=p, alternative='less')
>>> ci = res.confidence_interval(confidence_level=alpha)
>>> for x in rvs[rvs <= ci.high]:
...     res = stats.quantile_test(rvs, q=x, p=p, alternative='less')
...     assert res.pvalue > 1-alpha
>>> for x in rvs[rvs > ci.high]:
...     res = stats.quantile_test(rvs, q=x, p=p, alternative='less')
...     assert res.pvalue < 1-alpha 

此外,如果针对随机样本重复生成 95%置信区间,则在大约 95%的复制中,置信区间将包含真实的分位值。

>>> dist = stats.rayleigh() # our "unknown" distribution
>>> p = 0.2
>>> true_stat = dist.ppf(p) # the true value of the statistic
>>> n_trials = 1000
>>> quantile_ci_contains_true_stat = 0
>>> for i in range(n_trials):
...     data = dist.rvs(size=100, random_state=rng)
...     res = stats.quantile_test(data, p=p)
...     ci = res.confidence_interval(0.95)
...     if ci[0] < true_stat < ci[1]:
...         quantile_ci_contains_true_stat += 1
>>> quantile_ci_contains_true_stat >= 950
True 

只要样本是独立同分布的,这对任何分布和任何分位数都适用。

scipy.stats.skewtest

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

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

检验偏度是否与正态分布不同。

此函数测试空假设,即样本抽取的总体偏度与对应的正态分布的偏度相同。

参数:

aarray

待测试的数据。

axisint 或 None,可选

计算统计量的轴线。默认为 0。如果为 None,则在整个数组 a 上计算。

nan_policy,可选

定义当输入包含 NaN 时如何处理。有以下选项可用(默认为 'propagate'):

  • ‘propagate’:返回 NaN

  • ‘raise’:抛出错误

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

alternative,可选

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

  • ‘two-sided’:样本背后的分布的偏度与正态分布(即 0)不同

  • ‘less’:样本背后的分布的偏度小于正态分布的偏度

  • ‘greater’:样本背后的分布的偏度大于正态分布的偏度

1.7.0 版中的新功能。

返回:

statisticfloat

此测试的计算 z 分数。

pvaluefloat

假设检验的 p 值。

注意

样本大小必须至少为 8。

参考文献

[1]

R. B. D’Agostino, A. J. Belanger 和 R. B. D’Agostino Jr.,“使用强大和信息丰富的正态性检验的建议”,《美国统计学家》44,第 316-321 页,1990 年。

[2]

Shapiro,S. S.,& Wilk,M. B.(1965)。正态性的方差分析检验(完整样本)。《生物统计学》,52(3/4),591-611。

[3]

B. Phipson 和 G. K. Smyth。“当排列随机抽取时,排列 p 值永远不应为零:计算精确 p 值。”《遗传学和分子生物学中的统计应用》9.1(2010)。

示例

假设我们希望从测量中推断出,医学研究中成年男性的体重是否不符合正态分布 [2]。以下是记录在数组 x 中的体重(磅)。

>>> import numpy as np
>>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236]) 

来自[1]的偏度检验从计算基于样本偏度的统计量开始。

>>> from scipy import stats
>>> res = stats.skewtest(x)
>>> res.statistic
2.7788579769903414 

因为正态分布的偏度为零,所以此统计量的大小 tend 对于从正态分布中抽取的样本而言通常较低。

通过比较观察到的统计量的值与空假设的空分布来执行测试:统计量值的分布是在空假设下得出的,即权重是从正态分布中抽取的分布。

对于此测试,非常大样本的统计量的空分布是标准正态分布。

>>> import matplotlib.pyplot as plt
>>> dist = stats.norm()
>>> st_val = np.linspace(-5, 5, 100)
>>> pdf = dist.pdf(st_val)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def st_plot(ax):  # we'll reuse this
...     ax.plot(st_val, pdf)
...     ax.set_title("Skew Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> st_plot(ax)
>>> plt.show() 

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

比较由 p 值量化:在零分布中比观察到的统计值更极端或更极端的值的比例。在双侧检验中,零分布中大于观察统计量的元素和小于观察统计量的负值都被认为是“更极端”的。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> st_plot(ax)
>>> pvalue = dist.cdf(-res.statistic) + 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, (3, 0.005), (3.25, 0.02), arrowprops=props)
>>> i = st_val >= res.statistic
>>> ax.fill_between(st_val[i], y1=0, y2=pdf[i], color='C0')
>>> i = st_val <= -res.statistic
>>> ax.fill_between(st_val[i], y1=0, y2=pdf[i], color='C0')
>>> ax.set_xlim(-5, 5)
>>> ax.set_ylim(0, 0.1)
>>> plt.show() 

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

>>> res.pvalue
0.005455036974740185 

如果 p 值“小” - 即从一个正态分布的总体中抽取这样一个统计极值的概率很低 - 这可能被视为反对零假设的证据,赞成备择假设:这些权重不是从正态分布中抽取的。请注意:

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

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

注意标准正态分布提供了零分布的渐近逼近;对于有许多观测样本的情况,它只是准确的。对于像我们这样的小样本,scipy.stats.monte_carlo_test可能提供了一个更准确的,尽管是随机的,精确 p 值的近似。

>>> def statistic(x, axis):
...     # get just the skewtest statistic; ignore the p-value
...     return stats.skewtest(x, axis=axis).statistic
>>> res = stats.monte_carlo_test(x, stats.norm.rvs, statistic)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> st_plot(ax)
>>> ax.hist(res.null_distribution, np.linspace(-5, 5, 50),
...         density=True)
>>> ax.legend(['aymptotic approximation\n(many observations)',
...            'Monte Carlo approximation\n(11 observations)'])
>>> plt.show() 

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

>>> res.pvalue
0.0062  # may vary 

在这种情况下,渐近逼近和蒙特卡罗逼近非常接近,即使对于我们的小样本也是如此。

scipy.stats.kurtosistest

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

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

测试数据集是否具有正态峰度。

此函数检验假设,样本抽取自的总体峰度是正态分布的峰度。

参数:

a数组

样本数据的数组。

axis整数或 None, 可选

计算检验的轴线。默认为 0。如果为 None,则计算整个数组 a

nan_policy, 可选

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

  • ‘propagate’:返回 nan

  • ‘raise’:抛出错误

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

alternative, 可选

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

  • ‘two-sided’: 样本所在分布的峰度与正态分布不同。

  • ‘less’:样本所在分布的峰度小于正态分布的峰度

  • ‘greater’:样本所在分布的峰度大于正态分布的峰度

从版本 1.7.0 开始。

返回:

statistic浮点数

此检验的计算 z 分数。

pvalue浮点数

假设检验的 p 值。

注释

仅对 n>20 有效。此函数使用 [1] 中描述的方法。

参考文献

[1] (1,2)

参见例如 F. J. Anscombe, W. J. Glynn,“正态样本 b2 峰度统计量的分布”, Biometrika, vol. 70, pp. 227-234, 1983。

[2]

Shapiro, S. S., & Wilk, M. B. (1965). 方差分析检验正态性(完整样本)。 Biometrika, 52(3/4), 591-611.

[3]

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

[4]

Panagiotakos, D. B. (2008). 生物医学研究中 p 值的价值. The open cardiovascular medicine journal, 2, 97.

示例

假设我们希望从测量中推断,医学研究中成年男性的体重不服从正态分布[2]。以下是记录在数组 x 中的体重(磅)。

>>> import numpy as np
>>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236]) 

[1,2] 中的峰度检验首先基于样本(过量/费舍尔)峰度计算统计量。

>>> from scipy import stats
>>> res = stats.kurtosistest(x)
>>> res.statistic
2.3048235214240873 

(该检验警告我们的样本观察值太少,无法进行检验。我们将在示例结束时返回这一点。) 因为正态分布的过量峰度为零(定义如此),所以从正态分布中抽取的样本的此统计量的大小趋于较低。

测试是通过比较统计量的观察值与零分布进行的:在零假设下,权重是从正态分布中抽取的统计量值的分布。

对于这个测试,对于非常大的样本,统计量的零分布是标准正态分布。

>>> import matplotlib.pyplot as plt
>>> dist = stats.norm()
>>> kt_val = np.linspace(-5, 5, 100)
>>> pdf = dist.pdf(kt_val)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def kt_plot(ax):  # we'll reuse this
...     ax.plot(kt_val, pdf)
...     ax.set_title("Kurtosis Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> kt_plot(ax)
>>> plt.show() 

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

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

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> kt_plot(ax)
>>> pvalue = dist.cdf(-res.statistic) + 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, (3, 0.005), (3.25, 0.02), arrowprops=props)
>>> i = kt_val >= res.statistic
>>> ax.fill_between(kt_val[i], y1=0, y2=pdf[i], color='C0')
>>> i = kt_val <= -res.statistic
>>> ax.fill_between(kt_val[i], y1=0, y2=pdf[i], color='C0')
>>> ax.set_xlim(-5, 5)
>>> ax.set_ylim(0, 0.1)
>>> plt.show() 

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

>>> res.pvalue
0.0211764592113868 

如果 p 值“小” - 也就是说,从一个正态分布的人群中取样数据,产生如此极端的统计值的概率很低 - 这可能被视为反对零假设的证据,支持备择假设:权重并非来自正态分布。请注意:

  • 反之则不成立;也就是说,该检验不能用来支持零假设。

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

请注意,标准正态分布提供了零分布的渐近逼近;它仅适用于具有许多观察值的样本。这就是为什么我们在示例开始时收到了警告的原因;我们的样本非常小。在这种情况下,scipy.stats.monte_carlo_test可能提供更精确的,尽管是随机的精确 p 值的近似。

>>> def statistic(x, axis):
...     # get just the skewtest statistic; ignore the p-value
...     return stats.kurtosistest(x, axis=axis).statistic
>>> res = stats.monte_carlo_test(x, stats.norm.rvs, statistic)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> kt_plot(ax)
>>> ax.hist(res.null_distribution, np.linspace(-5, 5, 50),
...         density=True)
>>> ax.legend(['aymptotic approximation\n(many observations)',
...            'Monte Carlo approximation\n(11 observations)'])
>>> plt.show() 

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

>>> res.pvalue
0.0272  # may vary 

此外,尽管它们具有随机性质,以这种方式计算的 p 值可以用来精确控制零假设的错误拒绝率 [4]

scipy.stats.normaltest

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

scipy.stats.normaltest(a, axis=0, nan_policy='propagate')

测试样本是否不同于正态分布。

此函数测试样本是否来自正态分布的零假设。基于 D’Agostino 和 Pearson 的[1], [2] 检验,结合偏度和峰度产生正态性的全能检验。

参数:

aarray_like

包含待测试样本的数组。

axisint 或 None,可选

计算测试的轴。默认为 0。如果为 None,则在整个数组 a 上计算。

nan_policy,可选

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

  • ‘propagate’: 返回 nan
  • ‘raise’: 抛出错误
  • ‘omit’: 在计算中忽略 nan 值

返回:

statisticfloat 或 数组

s² + k²,其中 s 是由 skewtest 返回的 z 分数,k 是由 kurtosistest 返回的 z 分数。

pvaluefloat 或 数组

双侧卡方检验的 p 值。

参考文献

[1] (1,2)

D’Agostino, R. B. (1971), “适用于中等和大样本量的正态性全能检验”, 生物统计学, 58, 341-348

[2] (1,2)

D’Agostino, R. 和 Pearson, E. S. (1973), “偏离正态性的检验”, 生物统计学, 60, 613-622

[3]

Shapiro, S. S., & Wilk, M. B. (1965). 完整样本的方差分析正态性检验。生物统计学, 52(3/4), 591-611.

[4]

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

[5]

Panagiotakos, D. B. (2008). 生物医学研究中 p 值的价值。开放心血管医学期刊, 2, 97.

示例

假设我们希望根据测量推断成年人类男性体重在医学研究中是否不服从正态分布[3]。下面的数组 x 记录了体重(磅)。

>>> import numpy as np
>>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236]) 

[1][2] 的正态性检验首先基于样本偏度和峰度计算统计量。

>>> from scipy import stats
>>> res = stats.normaltest(x)
>>> res.statistic
13.034263121192582 

(测试警告我们的样本观测量太少,无法进行测试。我们将在示例结束时回到这个问题。) 因为正态分布具有零偏度和零(“过剩”或“费舍尔”)峰度,所以对于从正态分布中抽取的样本,该统计量的值趋向较低。

测试通过将统计量的观察值与空分布进行比较来执行:在空假设下,重量来自正态分布的统计值分布。对于这种正态性检验,对于非常大的样本,空分布是自由度为两的卡方分布。

>>> import matplotlib.pyplot as plt
>>> dist = stats.chi2(df=2)
>>> stat_vals = np.linspace(0, 16, 100)
>>> pdf = dist.pdf(stat_vals)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def plot(ax):  # we'll reuse this
...     ax.plot(stat_vals, pdf)
...     ax.set_title("Normality Test Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> plot(ax)
>>> plt.show() 

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

比较由 p 值量化:空分布中大于或等于观察到的统计量值的比例。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> pvalue = dist.sf(res.statistic)
>>> annotation = (f'p-value={pvalue:.6f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (13.5, 5e-4), (14, 5e-3), arrowprops=props)
>>> i = stat_vals >= res.statistic  # index more extreme statistic values
>>> ax.fill_between(stat_vals[i], y1=0, y2=pdf[i])
>>> ax.set_xlim(8, 16)
>>> ax.set_ylim(0, 0.01)
>>> plt.show() 

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

>>> res.pvalue
0.0014779023013100172 

如果 p 值“较小” - 即从一个正态分布的总体中抽取数据得到该统计量的极端值的概率较低 - 这可能被视为支持备择假设而非原假设的证据:重量并非来自正态分布。注意:

  • 反之不成立;即该检验不能用于提供支持空假设的证据。

  • 将被视为“小”的值的阈值是在分析数据之前应该做出的选择 [4],考虑到误拒空假设(错误地拒绝空假设)和假阴性(未能拒绝错误的空假设)的风险。

请注意,卡方分布提供空分布的渐近近似;它仅对观测值较多的样本准确。这就是我们在示例开始时收到警告的原因;我们的样本相当小。在这种情况下,scipy.stats.monte_carlo_test 可能提供更准确的、虽然是随机的、对精确 p 值的近似。

>>> def statistic(x, axis):
...     # Get only the `normaltest` statistic; ignore approximate p-value
...     return stats.normaltest(x, axis=axis).statistic
>>> res = stats.monte_carlo_test(x, stats.norm.rvs, statistic,
...                              alternative='greater')
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> ax.hist(res.null_distribution, np.linspace(0, 25, 50),
...         density=True)
>>> ax.legend(['aymptotic approximation (many observations)',
...            'Monte Carlo approximation (11 observations)'])
>>> ax.set_xlim(0, 14)
>>> plt.show() 

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

>>> res.pvalue
0.0082  # may vary 

此外,尽管其随机性,以这种方式计算的 p 值可以用来精确控制空假设的误拒率 [5]

scipy.stats.jarque_bera

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

scipy.stats.jarque_bera(x, *, axis=None, nan_policy='propagate', keepdims=False)

在样本数据上执行 Jarque-Bera 拟合度检验。

Jarque-Bera 检验测试样本数据的偏度和峰度是否与正态分布匹配。

请注意,此测试仅适用于足够大量的数据样本(>2000),因为测试统计量在渐近上具有自由度为 2 的卡方分布。

参数:

x类似数组

随机变量的观察。

整数或无,默认:无

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

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdims布尔值,默认:False

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

返回:

结果显著性结果

具有以下属性的对象:

统计量浮点数

测试统计量。

p 值浮点数

假设检验的 p 值。

注解

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

参考文献

[1]

Jarque, C. 和 Bera, A. (1980) “回归残差的正态性、等方差性和序列独立性的有效检验”,6 Econometric Letters 255-259.

[2]

Shapiro, S. S., & Wilk, M. B. (1965). 正态性的方差分析检验(完整样本)。Biometrika,52(3/4),591-611。

[3]

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

[4]

Panagiotakos, D. B. (2008). 生物医学研究中 p 值的价值。The open cardiovascular medicine journal, 2, 97.

示例

假设我们希望从测量中推断,医学研究中成年男性体重是否不服从正态分布[2]。下面数组x记录了体重(磅)。

>>> import numpy as np
>>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236]) 

Jarque-Bera 测试首先计算基于样本偏度和峰度的统计量。

>>> from scipy import stats
>>> res = stats.jarque_bera(x)
>>> res.statistic
6.982848237344646 

因为正态分布具有零偏度和零(“过量”或“Fisher”)峰度,这个统计量的值对于从正态分布中抽取的样本 tend to 较低。

测试通过比较统计量的观察值与零分布来执行:该统计值派生的零假设下的统计值分布,即权重来自正态分布。对于 Jarque-Bera 测试,非常大样本的零分布是具有两自由度的卡方分布。

>>> import matplotlib.pyplot as plt
>>> dist = stats.chi2(df=2)
>>> jb_val = np.linspace(0, 11, 100)
>>> pdf = dist.pdf(jb_val)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def jb_plot(ax):  # we'll reuse this
...     ax.plot(jb_val, pdf)
...     ax.set_title("Jarque-Bera Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> jb_plot(ax)
>>> plt.show() 

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

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

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> jb_plot(ax)
>>> pvalue = dist.sf(res.statistic)
>>> annotation = (f'p-value={pvalue:.6f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (7.5, 0.01), (8, 0.05), arrowprops=props)
>>> i = jb_val >= res.statistic  # indices of more extreme statistic values
>>> ax.fill_between(jb_val[i], y1=0, y2=pdf[i])
>>> ax.set_xlim(0, 11)
>>> ax.set_ylim(0, 0.3)
>>> plt.show() 

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

>>> res.pvalue
0.03045746622458189 

如果 p 值“小” - 也就是说,从正态分布人群中抽样产生统计量极端值的概率较低 - 这可能被视为反对零假设的证据,支持备选假设:体重未来自正态分布。请注意:

  • 相反则不成立;即,该测试不用于提供支持零假设的证据。

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

注意,卡方分布提供了零分布的渐近逼近;它仅对有大量观测样本的样本准确。对于像我们这样的小样本,scipy.stats.monte_carlo_test可能提供了一个更准确的,尽管是随机的,准确的 p 值近似。

>>> def statistic(x, axis):
...     # underlying calculation of the Jarque Bera statistic
...     s = stats.skew(x, axis=axis)
...     k = stats.kurtosis(x, axis=axis)
...     return x.shape[axis]/6 * (s**2 + k**2/4)
>>> res = stats.monte_carlo_test(x, stats.norm.rvs, statistic,
...                              alternative='greater')
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> jb_plot(ax)
>>> ax.hist(res.null_distribution, np.linspace(0, 10, 50),
...         density=True)
>>> ax.legend(['aymptotic approximation (many observations)',
...            'Monte Carlo approximation (11 observations)'])
>>> plt.show() 

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

>>> res.pvalue
0.0097  # may vary 

此外,尽管它们具有随机性质,通过这种方式计算的 p 值可以用来精确控制零假设的虚假拒绝率[4]

scipy.stats.shapiro

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

scipy.stats.shapiro(x)

执行 Shapiro-Wilk 正态性检验。

Shapiro-Wilk 测试检验的零假设是数据来自正态分布。

参数:

xarray_like

样本数据数组。

返回:

统计量float

检验统计量。

P 值float

假设检验的 P 值。

另请参见

anderson

Anderson-Darling 正态性检验

kstest

Kolmogorov-Smirnov 拟合优度检验。

注意事项

描述的算法见[4],但未实现所述的参数审查。对于 N > 5000,W 测试统计量是准确的,但 P 值可能不准确。

参考文献

[1]

www.itl.nist.gov/div898/handbook/prc/section2/prc213.htm DOI:10.18434/M32189

[2] (1,2)

Shapiro, S. S. & Wilk, M.B, “方差分析正态性检验(完全样本)”,1965 年,《生物统计学》,第 52 卷,第 591-611 页,DOI:10.2307/2333709

[3]

Razali, N. M. & Wah, Y. B., “Shapiro-Wilk、Kolmogorov-Smirnov、Lilliefors 和 Anderson-Darling 测试的功效比较”,2011 年,《统计建模与分析杂志》,第 2 卷,第 21-33 页。

[4]

Royston P., “备注 AS R94: 关于算法 AS 181 的备注:W-测试用于正态性”,1995 年,《应用统计学》,第 44 卷,DOI:10.2307/2986146

[5]

Phipson B., and Smyth, G. K., “当排列随机抽取时,置换 P 值不应为零:计算精确 P 值”,2010 年,《统计应用于遗传学和分子生物学》,第 9 卷,DOI:10.2202/1544-6115.1585

[6]

Panagiotakos, D. B., “在生物医学研究中 P 值的价值”,《开放心血管医学杂志》,2008 年,第 2 卷,第 97-99 页,DOI:10.2174/1874192400802010097

示例

假设我们希望根据测量数据推断医学研究中成年男性体重是否不服从正态分布[2]。以下是记录在数组x中的成年男性体重(磅)。

>>> import numpy as np
>>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236]) 

[1][2]的正态性检验首先基于观测值与正态分布的预期顺序统计量之间的关系计算统计量。

>>> from scipy import stats
>>> res = stats.shapiro(x)
>>> res.statistic
0.7888147830963135 

该统计量的值对于从正态分布中抽取的样本趋于高(接近 1)。

该检验通过比较统计量的观察值与零分布进行:统计值形成的分布在零假设下,即权重是从正态分布中抽取的。对于这个正态性检验,零分布不容易精确计算,因此通常通过蒙特卡罗方法来近似,即从一个与x相同大小的正态分布中抽取许多样本,并计算每个统计量的值。

>>> def statistic(x):
...     # Get only the `shapiro` statistic; ignore its p-value
...     return stats.shapiro(x).statistic
>>> ref = stats.monte_carlo_test(x, stats.norm.rvs, statistic,
...                              alternative='less')
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> bins = np.linspace(0.65, 1, 50)
>>> def plot(ax):  # we'll reuse this
...     ax.hist(ref.null_distribution, density=True, bins=bins)
...     ax.set_title("Shapiro-Wilk Test Null Distribution \n"
...                  "(Monte Carlo Approximation, 11 Observations)")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> plot(ax)
>>> plt.show() 

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

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

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> plot(ax)
>>> annotation = (f'p-value={res.pvalue:.6f}\n(highlighted area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (0.75, 0.1), (0.68, 0.7), arrowprops=props)
>>> i_extreme = np.where(bins <= res.statistic)[0]
>>> for i in i_extreme:
...     ax.patches[i].set_color('C1')
>>> plt.xlim(0.65, 0.9)
>>> plt.ylim(0, 4)
>>> plt.show
>>> res.pvalue
0.006703833118081093 

如果 p 值“小” - 即从一个正态分布的总体中抽样数据,使得统计量的极端值的概率很低 - 这可能被视为反对零假设,支持备选假设:权重并非来自正态分布。注意:

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

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

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

scipy.stats.anderson

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

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

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

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

参数:

xarray_like

样本数据数组。

dist, 可选

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

返回:

resultAndersonResult

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

statisticfloat

Anderson-Darling 检验统计量。

critical_valueslist

此分布的关键值。

significance_levellist

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

fit_resultFitResult

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

另请参见

kstest

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

注释

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

正态/指数

15%,10%,5%,2.5%,1%

逻辑分布

25%,10%,5%,2.5%,1%,0.5%

gumbel_l / gumbel_r

25%,10%,5%,2.5%,1%

威布尔最小型

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

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

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

参考文献

[1]

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

[2]

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

[3]

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

[4]

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

[5]

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

[6]

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

[7]

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

例子

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

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

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

scipy.stats.cramervonmises

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

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

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

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

参数:

rvsarray_like

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

cdfstr 或 可调用对象

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

args元组,可选

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

axis整数或 None,默认为 0

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

nan_policy

定义如何处理输入的 NaN。

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

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

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

keepdims布尔值,默认为 False

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

返回:

res具有属性的对象

统计量为 float

Cramér-von Mises 统计量。

p 值为 float

p 值。

参见

kstest, cramervonmises_2samp

注释

从版本 1.6.0 开始。

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

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

参考文献

[1]

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

[2]

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

示例

假设我们希望测试由 scipy.stats.norm.rvs 生成的数据是否实际上是从标准正态分布中抽取的。我们选择显著性水平 alpha=0.05

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> x = stats.norm.rvs(size=500, random_state=rng)
>>> res = stats.cramervonmises(x, 'norm')
>>> res.statistic, res.pvalue
(0.1072085112565724, 0.5508482238203407) 

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

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

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

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

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

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

scipy.stats.ks_1samp

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

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

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

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

参数:

x array_like

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

cdf 可调用函数

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

args 元组,序列,可选

cdf 一起使用的分布参数。

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

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

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

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

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

axis 整数或 None,默认为 0

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

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

定义如何处理输入的 NaN。

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

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

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

keepdims 布尔值,默认为 False

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

返回:

res:KstestResult

一个包含属性的对象:

statisticfloat

KS 检验统计量,可以是 D+、D- 或 D(两者中的最大值)

pvaluefloat

单尾或双尾 p 值。

statistic_locationfloat

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

statistic_signint

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

另请参见

ks_2sampkstest

注意

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

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

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

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

注意备择假设描述的是底层分布的CDF,而不是观察值。例如,假设 x1 ~ F 和 x2 ~ G。如果对所有 x,F(x) > G(x),那么 x1 中的值往往小于 x2 中的值。

从 SciPy 1.9 开始,不推荐使用np.matrix输入进行计算前会被转换为np.ndarray。在这种情况下,输出将是一个标量或适当形状的np.ndarray,而不是 2D 的np.matrix。类似地,忽略掩码数组的掩码元素时,输出将是一个标量或np.ndarray,而不是带有mask=False的掩码数组。

例子

假设我们希望测试一个样本是否符合标准正态分布的零假设。我们选择 95%的置信水平;也就是说,如果 p 值小于 0.05,我们将拒绝零假设,支持备择假设。

在测试均匀分布数据时,我们预期将会拒绝零假设。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> stats.ks_1samp(stats.uniform.rvs(size=100, random_state=rng),
...                stats.norm.cdf)
KstestResult(statistic=0.5001899973268688, pvalue=1.1616392184763533e-23) 

实际上,p 值低于我们的 0.05 阈值,因此我们拒绝零假设,支持默认的“双侧”备择假设:数据并按照标准正态分布分布。

在测试标准正态分布的随机变量时,我们期望数据大部分时间与零假设一致。

>>> x = stats.norm.rvs(size=100, random_state=rng)
>>> stats.ks_1samp(x, stats.norm.cdf)
KstestResult(statistic=0.05345882212970396, pvalue=0.9227159037744717) 

正如预期的那样,p 值为 0.92 不低于我们的 0.05 阈值,因此我们无法拒绝零假设。

然而,假设随机变量分布于一个向更大数值偏移的正态分布。在这种情况下,底层分布的累积密度函数(CDF)倾向于比标准正态分布的 CDF更少。因此,我们预期会以alternative='less'的方式拒绝零假设:

>>> x = stats.norm.rvs(size=100, loc=0.5, random_state=rng)
>>> stats.ks_1samp(x, stats.norm.cdf, alternative='less')
KstestResult(statistic=0.17482387821055168, pvalue=0.001913921057766743) 

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

scipy.stats.goodness_of_fit

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

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

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

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

参数:

distscipy.stats.rv_continuous

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

data1D array_like

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

known_paramsdict,可选

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

fit_paramsdict,可选

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

guessed_paramsdict,可选

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

statistic,可选

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

n_mc_samplesint,默认值:9999

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

random_state{None, int, numpy.random.Generator

numpy.random.RandomState,可选

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

如果random_stateNone(默认),则使用numpy.random.RandomState单例。如果random_state为整数,则使用新的RandomState实例,并使用random_state作为种子。如果random_state已经是GeneratorRandomState实例,则使用提供的实例。

返回:

resGoodnessOfFitResult

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

fit_resultFitResult

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

statisticfloat

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

pvaluefloat

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

null_distributionndarray

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

注意事项

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

传统的拟合优度检验

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

算法概述

简言之,此例程执行以下步骤:

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

更详细地说,步骤如下。

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

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

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

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

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

限制

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

反模式

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

参考文献

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

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

[2]

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

[3]

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

[4]

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

[5]

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

[6]

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

[7]

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

示例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

这个图看起来令人放心。

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

>>> res.statistic, res.pvalue
(0.2231991510248692, 0.0525) 
posted @ 2024-06-27 17:06  绝不原创的飞龙  阅读(2)  评论(0编辑  收藏  举报