SciPy 1.12 中文文档（六十六）
scipy.stats.boxcox_normmax
scipy.stats.boxcox_normmax(x, brack=None, method='pearsonr', optimizer=None)
计算输入数据的最佳 Box-Cox 变换参数。
参数：
x array_like
输入数组。所有条目必须为正有限实数。
brack 2-元组，可选，默认为 (-2.0, 2.0)
默认 optimize.brent 求解器进行向下斜坡搜索的起始区间。请注意，这在大多数情况下并不重要；最终结果允许在此区间之外。如果传递了 optimizer，则 brack 必须为 None。
method str，可选
确定最佳变换参数（boxcox 的 lmbda 参数）的方法。选项包括：
‘pearsonr’（默认）
最大化皮尔逊相关系数，使 y = boxcox(x) 与 y 的期望值在 x 若为正态分布时保持一致。
‘mle’
最大化对数似然 boxcox_llf。这是 boxcox 中使用的方法。
‘all’
使用所有可用的优化方法，并返回所有结果。有助于比较不同的方法。
optimizer 可调用函数，可选
optimizer 是一个接受一个参数的可调用函数：
funcallable
要最小化的目标函数。fun 接受一个参数，即 Box-Cox 变换参数 lmbda，并返回在提供的参数处函数值（例如，负对数似然）。optimizer 的任务是找到最小化 fun 的 lmbda 值。
并返回一个对象，如 scipy.optimize.OptimizeResult 的实例，其中保存了属性 x 的最佳 lmbda 值。
更多信息请参见下面的示例或 scipy.optimize.minimize_scalar 的文档。
返回：
maxlog 浮点数或 ndarray
找到的最佳变换参数。对于 method='all'，这是一个数组而不是标量。
参见
boxcox、boxcox_llf、boxcox_normplot、scipy.optimize.minimize_scalar 的示例
示例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
我们可以生成一些数据，并以各种方式确定最佳的 lmbda：
>>> rng = np.random.default_rng()
>>> x = stats.loggamma.rvs(5, size=30, random_state=rng) + 5
>>> y, lmax_mle = stats.boxcox(x)
>>> lmax_pearsonr = stats.boxcox_normmax(x)
>>> lmax_mle
2.217563431465757
>>> lmax_pearsonr
2.238318660200961
>>> stats.boxcox_normmax(x, method='all')
array([2.23831866, 2.21756343])
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> prob = stats.boxcox_normplot(x, -10, 10, plot=ax)
>>> ax.axvline(lmax_mle, color='r')
>>> ax.axvline(lmax_pearsonr, color='g', ls='--')
>>> plt.show()
或者，我们可以定义自己的优化器函数。假设我们只对区间[6, 7]上的lmbda值感兴趣，我们希望使用scipy.optimize.minimize_scalar，设置method='bounded'，并且在优化对数似然函数时希望使用更严格的公差。为此，我们定义一个接受位置参数fun的函数，并使用scipy.optimize.minimize_scalar来最小化fun，同时满足提供的边界和公差：
>>> from scipy import optimize
>>> options = {'xatol': 1e-12} # absolute tolerance on `x`
>>> def optimizer(fun):
... return optimize.minimize_scalar(fun, bounds=(6, 7),
... method="bounded", options=options)
>>> stats.boxcox_normmax(x, optimizer=optimizer)
6.000...
scipy.stats.boxcox_llf
scipy.stats.boxcox_llf(lmb, data)
Box-Cox 对数似然函数。
参数：
lmb 标量
Box-Cox 转换的参数。详情请见 boxcox。
data array_like
计算 Box-Cox 对数似然函数的数据。如果 data 是多维的，则沿着第一轴计算对数似然。
返回：
llf 浮点数或者 ndarray
给定 lmb 的 Box-Cox 对数似然函数。对于 1-D data 是一个浮点数，否则是一个数组。
参见
boxcox, probplot, boxcox_normplot, boxcox_normmax
注意事项
Box-Cox 对数似然函数在这里定义为
[llf = (\lambda - 1) \sum_i(\log(x_i)) - N/2 \log(\sum_i (y_i - \bar{y})² / N),]
其中 y 是经过 Box-Cox 变换的输入数据 x。
示例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
生成一些随机变量并计算它们的 Box-Cox 对数似然值，使用一系列 lmbda 值：
>>> rng = np.random.default_rng()
>>> x = stats.loggamma.rvs(5, loc=10, size=1000, random_state=rng)
>>> lmbdas = np.linspace(-2, 10)
>>> llf = np.zeros(lmbdas.shape, dtype=float)
>>> for ii, lmbda in enumerate(lmbdas):
... llf[ii] = stats.boxcox_llf(lmbda, x)
同时使用 boxcox 找到最优 lmbda 值：
>>> x_most_normal, lmbda_optimal = stats.boxcox(x)
绘制以 lmbda 为函数的对数似然函数图。添加最优 lmbda 作为水平线来确认确实是最优值：
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(lmbdas, llf, 'b.-')
>>> ax.axhline(stats.boxcox_llf(lmbda_optimal, x), color='r')
>>> ax.set_xlabel('lmbda parameter')
>>> ax.set_ylabel('Box-Cox log-likelihood')
现在添加一些概率图，以展示通过 boxcox 变换的数据在对数似然函数最大化时看起来最接近正态分布：
>>> locs = [3, 10, 4] # 'lower left', 'center', 'lower right'
>>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
... xt = stats.boxcox(x, lmbda=lmbda)
... (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
... ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
... ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
... ax_inset.set_xticklabels([])
... ax_inset.set_yticklabels([])
... ax_inset.set_title(r'$\lambda=%1.2f$' % lmbda)
>>> plt.show()
scipy.stats.yeojohnson
scipy.stats.yeojohnson(x, lmbda=None)
返回经 Yeo-Johnson 功率变换后的数据集。
参数：
xndarray
输入数组。应为一维数组。
lmbdafloat，可选
如果 lmbda 为 None，则找到最大化对数似然函数的 lambda，并将其作为第二个输出参数返回。否则，按给定值进行变换。
返回：
yeojohnson：ndarray
经 Yeo-Johnson 功率变换后的数组。
maxlogfloat，可选
如果 lmbda 参数为 None，则第二个返回参数为最大化对数似然函数的 lambda。
参见
probplot，yeojohnson_normplot，yeojohnson_normmax，yeojohnson_llf，boxcox
注：
Yeo-Johnson 变换由以下式给出：
y = ((x + 1)**lmbda - 1) / lmbda, for x >= 0, lmbda != 0
log(x + 1), for x >= 0, lmbda = 0
-((-x + 1)**(2 - lmbda) - 1) / (2 - lmbda), for x < 0, lmbda != 2
-log(-x + 1), for x < 0, lmbda = 2
与 boxcox 不同，yeojohnson 不要求输入数据为正数。
自 1.2.0 版新增。
参考文献
I. Yeo 和 R.A. Johnson，《改善正态性或对称性的新型功率变换家族》，Biometrika 87.4 (2000)：
示例
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
我们从非正态分布生成一些随机变量，并为其制作概率图，以显示其在尾部不是正态分布：
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(211)
>>> x = stats.loggamma.rvs(5, size=500) + 5
>>> prob = stats.probplot(x, dist=stats.norm, plot=ax1)
>>> ax1.set_xlabel('')
>>> ax1.set_title('Probplot against normal distribution')
我们现在使用 yeojohnson 对数据进行变换，使其最接近正态分布：
>>> ax2 = fig.add_subplot(212)
>>> xt, lmbda = stats.yeojohnson(x)
>>> prob = stats.probplot(xt, dist=stats.norm, plot=ax2)
>>> ax2.set_title('Probplot after Yeo-Johnson transformation')
>>> plt.show()
scipy.stats.yeojohnson_normmax
scipy.stats.yeojohnson_normmax(x, brack=None)
计算最优的 Yeo-Johnson 变换参数。
计算输入数据的最优 Yeo-Johnson 变换参数，使用最大似然估计。
参数：
x 类似数组
输入数组。
brack 2 元组，可选
用于 optimize.brent 的下坡搜索的起始区间。请注意，在大多数情况下这并不关键；最终结果允许超出此区间。如果为 None，则使用 optimize.fminbound 并设置避免溢出的边界。
返回值：
maxlog 浮点数
找到的最优变换参数。
另请参阅
yeojohnson, yeojohnson_llf, yeojohnson_normplot
注意事项
自版本 1.2.0 起新增。
示例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
生成一些数据并确定最优的 lmbda
>>> rng = np.random.default_rng()
>>> x = stats.loggamma.rvs(5, size=30, random_state=rng) + 5
>>> lmax = stats.yeojohnson_normmax(x)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> prob = stats.yeojohnson_normplot(x, -10, 10, plot=ax)
>>> ax.axvline(lmax, color='r')
>>> plt.show()
scipy.stats.yeojohnson_llf
scipy.stats.yeojohnson_llf(lmb, data)
Yeo-Johnson 对数似然函数。
参数：
lmb标量
Yeo-Johnson 变换的参数。详情请参阅yeojohnson。
数据array_like
用于计算 Yeo-Johnson 对数似然的数据。如果data是多维的，则沿第一轴计算对数似然。
返回：
llf浮点数
给定lmb的 Yeo-Johnson 对数似然函数。
另请参阅
yeojohnson，probplot，yeojohnson_normplot，yeojohnson_normmax
注意事项
Yeo-Johnson 对数似然函数在这里定义为
[llf = -N/2 \log(\hat{\sigma}²) + (\lambda - 1) \sum_i \text{ sign }(x_i)\log(|x_i| + 1)]
其中(\hat{\sigma}²)是 Yeo-Johnson 变换后输入数据x的估计方差。
版本 1.2.0 中的新功能。
举例
>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
生成一些随机变量，并计算它们的 Yeo-Johnson 对数似然值，用一系列lmbda值：
>>> x = stats.loggamma.rvs(5, loc=10, size=1000)
>>> lmbdas = np.linspace(-2, 10)
>>> llf = np.zeros(lmbdas.shape, dtype=float)
>>> for ii, lmbda in enumerate(lmbdas):
... llf[ii] = stats.yeojohnson_llf(lmbda, x)
还可以使用yeojohnson找到最优的lmbda值：
>>> x_most_normal, lmbda_optimal = stats.yeojohnson(x)
绘制对数似然函数作为lmbda的函数。添加最优lmbda作为水平线以检查是否确实是最优：
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(lmbdas, llf, 'b.-')
>>> ax.axhline(stats.yeojohnson_llf(lmbda_optimal, x), color='r')
>>> ax.set_xlabel('lmbda parameter')
>>> ax.set_ylabel('Yeo-Johnson log-likelihood')
现在添加一些概率图，显示对数似然函数最大化的地方，用yeojohnson变换后的数据看起来最接近正态分布：
>>> locs = [3, 10, 4] # 'lower left', 'center', 'lower right'
>>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
... xt = stats.yeojohnson(x, lmbda=lmbda)
... (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
... ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
... ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
... ax_inset.set_xticklabels([])
... ax_inset.set_yticklabels([])
... ax_inset.set_title(r'$\lambda=%1.2f$' % lmbda)
>>> plt.show()
scipy.stats.obrientransform
scipy.stats.obrientransform(*samples)
计算输入数据（任意数量的数组）上的 O’Brien 变换。
用于在运行单因素统计之前测试方差的均匀性。*samples中的每个数组都是因素的一个水平。如果在转换后的数据上运行 f_oneway，并且发现显著性，则方差不相等。来自 Maxwell 和 Delaney [1]，p.112。
参数：
**sample1, sample2, …**array_like
任意数量的数组。
返回：
obrientransformndarray
用于 ANOVA 的转换数据。结果的第一个维度对应于转换数组的序列。如果给定的数组都是相同长度的 1-D 数组，则返回值是一个 2-D 数组；否则它是一个对象类型的 1-D 数组，其中每个元素都是一个 ndarray。
参考文献
[1]
S. E. Maxwell 和 H. D. Delaney，“Designing Experiments and Analyzing Data: A Model Comparison Perspective”，Wadsworth，1990 年。
示例
我们将测试以下数据集的方差差异。
>>> x = [10, 11, 13, 9, 7, 12, 12, 9, 10]
>>> y = [13, 21, 5, 10, 8, 14, 10, 12, 7, 15]
对数据应用 O’Brien 变换。
>>> from scipy.stats import obrientransform
>>> tx, ty = obrientransform(x, y)
使用 scipy.stats.f_oneway 对转换数据应用单因素 ANOVA 检验。
>>> from scipy.stats import f_oneway
>>> F, p = f_oneway(tx, ty)
>>> p
0.1314139477040335
如果我们要求 p < 0.05 表示显著性，则我们不能断定方差不同。
scipy.stats.sigmaclip
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.sigmaclip.html#scipy.stats.sigmaclip
scipy.stats.sigmaclip(a, low=4.0, high=4.0)
对数组元素执行迭代的 Sigma 切除。
从完整样本开始，移除所有在临界范围之外的元素，即满足以下条件之一的输入数组 c 的所有元素：
c < mean(c) - std(c)*low
c > mean(c) + std(c)*high
迭代继续进行，直到没有元素在（更新后的）范围之外。
参数：
aarray_like
数据数组，如果不是 1-D，则会展平。
lowfloat，可选
Sigma 切除的下限系数。默认为 4。
highfloat，可选
Sigma 切除的上限系数。默认为 4。
返回：
clippedndarray
带有切除元素的输入数组。
lowerfloat
用于切除的下阈值。
upperfloat
用于切除的上阈值。
示例：
>>> import numpy as np
>>> from scipy.stats import sigmaclip
>>> a = np.concatenate((np.linspace(9.5, 10.5, 31),
... np.linspace(0, 20, 5)))
>>> fact = 1.5
>>> c, low, upp = sigmaclip(a, fact, fact)
>>> c
array([ 9.96666667, 10\. , 10.03333333, 10\. ])
>>> c.var(), c.std()
(0.00055555555555555165, 0.023570226039551501)
>>> low, c.mean() - fact*c.std(), c.min()
(9.9646446609406727, 9.9646446609406727, 9.9666666666666668)
>>> upp, c.mean() + fact*c.std(), c.max()
(10.035355339059327, 10.035355339059327, 10.033333333333333)
>>> a = np.concatenate((np.linspace(9.5, 10.5, 11),
... np.linspace(-100, -50, 3)))
>>> c, low, upp = sigmaclip(a, 1.8, 1.8)
>>> (c == np.linspace(9.5, 10.5, 11)).all()
True
scipy.stats.trimboth
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.trimboth.html#scipy.stats.trimboth
scipy.stats.trimboth(a, proportiontocut, axis=0)
从数组的两端切除一部分项目。
从传递的数组的两端切除传递的项目的比例（即，proportiontocut = 0.1，切片左边的 10% 和 右边的 10% 的分数）。修剪的值是最低和最高的值。如果比例导致非整数切片索引，则切片较少（即，保守地切片 proportiontocut）。
参数：
aarray_like
要修剪的数据。
proportiontocutfloat
要修剪的每端的总数据集的比例（范围在 0-1 之间）。
axisint 或 None，可选
数据修剪的轴。默认为 0。如果为 None，则在整个数组 a 上计算。
返回：
outndarray
数组 a 的修剪版本。修剪内容的顺序未定义。
另请参见
trim_mean
示例
创建一个包含 10 个值的数组，并从每端修剪 10% 的值：
>>> import numpy as np
>>> from scipy import stats
>>> a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> stats.trimboth(a, 0.1)
array([1, 3, 2, 4, 5, 6, 7, 8])
输入数组的元素根据值进行修剪，但输出数组未必按顺序排列。
要修剪的比例向下舍入到最接近的整数。例如，从一个包含 10 个值的数组的每端修剪 25% 的值将返回一个包含 6 个值的数组：
>>> b = np.arange(10)
>>> stats.trimboth(b, 1/4).shape
(6,)
可以沿任何轴或整个数组修剪多维数组：
>>> c = [2, 4, 6, 8, 0, 1, 3, 5, 7, 9]
>>> d = np.array([a, b, c])
>>> stats.trimboth(d, 0.4, axis=0).shape
(1, 10)
>>> stats.trimboth(d, 0.4, axis=1).shape
(3, 2)
>>> stats.trimboth(d, 0.4, axis=None).shape
(6,)
scipy.stats.trim1
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.trim1.html#scipy.stats.trim1
scipy.stats.trim1(a, proportiontocut, tail='right', axis=0)
从传递的数组分布的一端切片掉一部分。
如果 proportiontocut = 0.1，则切除分数‘最左侧’或‘最右侧’的 10%得分。修剪较少，如果比例导致非整数切片索引（即保守地切掉 proportiontocut ）。
参数：
a 数组样式
输入数组。
proportiontocut 浮点数
分数是从分布的“左侧”或“右侧”截掉的。
tail {‘left’，‘right’}，可选
默认为‘right’。
axis 整数或无，可选
用于修剪数据的轴。默认为 0。如果为 None，则在整个数组 a 上计算。
返回：
trim1 数组
缩短版本的数组 a。修剪后内容的顺序未定义。
示例
创建包含 10 个值的数组并修剪其最低值的 20%：
>>> import numpy as np
>>> from scipy import stats
>>> a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> stats.trim1(a, 0.2, 'left')
array([2, 4, 3, 5, 6, 7, 8, 9])
请注意，输入数组的元素按值修剪，但输出数组未必排序。
要修剪的比例向下舍入到最接近的整数。例如，从包含 10 个值的数组中修剪 25%的值将返回包含 8 个值的数组：
>>> b = np.arange(10)
>>> stats.trim1(b, 1/4).shape
(8,)
多维数组可以沿任意轴或整个数组进行修剪：
>>> c = [2, 4, 6, 8, 0, 1, 3, 5, 7, 9]
>>> d = np.array([a, b, c])
>>> stats.trim1(d, 0.8, axis=0).shape
(1, 10)
>>> stats.trim1(d, 0.8, axis=1).shape
(3, 2)
>>> stats.trim1(d, 0.8, axis=None).shape
(6,)
scipy.stats.zmap
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.zmap.html#scipy.stats.zmap
scipy.stats.zmap(scores, compare, axis=0, ddof=0, nan_policy='propagate')
计算相对 z-scores。
返回一个 z-score 数组，即标准化为零均值和单位方差的分数，其中均值和方差是从比较数组计算得出的。
参数：
scores：array_like
用于计算 z-scores 的输入。
compare：array_like
用于计算归一化均值和标准差的输入；假设与scores具有相同的维度。
axis：int 或 None，可选
计算compare的均值和方差的轴。默认为 0。如果为 None，则在整个数组scores上计算。
ddof：int，可选
在标准差计算中的自由度校正。默认为 0。
nan_policy：{‘propagate’, ‘raise’, ‘omit’}，可选
定义如何处理compare中 NaN 的出现。‘propagate’返回 NaN，‘raise’引发异常，‘omit’执行计算时忽略 NaN 值。默认为‘propagate’。请注意，当值为‘omit’时，scores中的 NaN 也会传播到输出，但它们不会影响对非 NaN 值计算的 z-scores。
返回：
zscore：array_like
与scores相同形状的 Z-scores。
注意
此函数保留 ndarray 子类，并且还适用于矩阵和掩码数组（它使用asanyarray而不是asarray作为参数）。
示例
>>> from scipy.stats import zmap
>>> a = [0.5, 2.0, 2.5, 3]
>>> b = [0, 1, 2, 3, 4]
>>> zmap(a, b)
array([-1.06066017, 0\. , 0.35355339, 0.70710678])
scipy.stats.zscore
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.zscore.html#scipy.stats.zscore
scipy.stats.zscore(a, axis=0, ddof=0, nan_policy='propagate')
计算 z 分数。
计算样本中每个值相对于样本均值和标准差的 z 分数。
参数：
aarray_like
一个类似数组的对象，包含样本数据。
axisint 或 None，可选
操作的轴。默认为 0。如果为 None，则在整个数组a上计算。
ddofint，可选
在标准差计算中的自由度修正。默认为 0。
nan_policy{‘propagate’, ‘raise’, ‘omit’}，可选
定义输入包含 nan 时的处理方式。‘propagate’返回 nan，‘raise’抛出错误，‘omit’执行计算时忽略 nan 值。默认为‘propagate’。注意，当值为‘omit’时，输入中的 nan 也会传播到输出，但它们不会影响计算非 nan 值的 z 分数。
返回：
zscorearray_like
标准化后的 z 分数，按输入数组a的均值和标准差计算。
另请参阅
算术平均
算术标准差
scipy.stats.gzscore
几何标准分数
注释
此函数保留 ndarray 子类，并且还适用于矩阵和掩码数组（它使用asanyarray而不是asarray作为参数）。
参考文献
[1]
“标准分数”，维基百科，zh.wikipedia.org/wiki/%E6%A8%99%E6%BA%96%E5%88%86%E6%95%B8。
[2]
Huck, S. W., Cross, T. L., Clark, S. B，“克服关于 Z 分数的误解”，《教学统计学》，第 8 卷，第 38-40 页，1986 年
示例
>>> import numpy as np
>>> a = np.array([ 0.7972, 0.0767, 0.4383, 0.7866, 0.8091,
... 0.1954, 0.6307, 0.6599, 0.1065, 0.0508])
>>> from scipy import stats
>>> stats.zscore(a)
array([ 1.1273, -1.247 , -0.0552, 1.0923, 1.1664, -0.8559, 0.5786,
0.6748, -1.1488, -1.3324])
沿指定轴计算，使用 n-1 自由度(ddof=1)计算标准差：
>>> b = np.array([[ 0.3148, 0.0478, 0.6243, 0.4608],
... [ 0.7149, 0.0775, 0.6072, 0.9656],
... [ 0.6341, 0.1403, 0.9759, 0.4064],
... [ 0.5918, 0.6948, 0.904 , 0.3721],
... [ 0.0921, 0.2481, 0.1188, 0.1366]])
>>> stats.zscore(b, axis=1, ddof=1)
array([[-0.19264823, -1.28415119, 1.07259584, 0.40420358],
[ 0.33048416, -1.37380874, 0.04251374, 1.00081084],
[ 0.26796377, -1.12598418, 1.23283094, -0.37481053],
[-0.22095197, 0.24468594, 1.19042819, -1.21416216],
[-0.82780366, 1.4457416 , -0.43867764, -0.1792603 ]])
以*nan_policy='omit'*为例：
>>> x = np.array([[25.11, 30.10, np.nan, 32.02, 43.15],
... [14.95, 16.06, 121.25, 94.35, 29.81]])
>>> stats.zscore(x, axis=1, nan_policy='omit')
array([[-1.13490897, -0.37830299, nan, -0.08718406, 1.60039602],
[-0.91611681, -0.89090508, 1.4983032 , 0.88731639, -0.5785977 ]])
scipy.stats.gzscore
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.gzscore.html#scipy.stats.gzscore
scipy.stats.gzscore(a, *, axis=0, ddof=0, nan_policy='propagate')
计算几何标准分数。
计算样本中每个严格正值的几何 z 分数，相对于几何平均值和标准差。数学上，可以将几何 z 分数评估为：
gzscore = log(a/gmu) / log(gsigma)
其中gmu（或gsigma）是几何平均值（或标准差）。
参数：
a类似数组
样本数据。
axis整数或 None，可选
操作的轴。默认为 0。如果为 None，则在整个数组a上计算。
ddof整数，可选
在标准差计算中的自由度校正。默认为 0。
nan_policy{‘propagate’, ‘raise’, ‘omit’}，可选
定义输入包含 NaN 时的处理方式。‘propagate’返回 NaN，‘raise’引发错误，‘omit’执行计算时忽略 NaN 值。默认为‘propagate’。注意，当值为‘omit’时，输入中的 NaN 也会传播到输出，但它们不会影响对非 NaN 值计算的几何 z 分数。
返回：
gzscore类似数组
几何 z 分数，通过输入数组a的几何平均值和几何标准差标准化。
另请参见
gmean
几何平均值
gstd
几何标准差
zscore
标准分数
注意
此函数保留 ndarray 子类，并且也适用于矩阵和掩码数组（它使用asanyarray而不是asarray作为参数）。
新版本中的新功能 1.8。
参考文献
[1]
“几何标准分数”，维基百科，en.wikipedia.org/wiki/Geometric_standard_deviation#Geometric_standard_score。
示例
从对数正态分布中抽样：
>>> import numpy as np
>>> from scipy.stats import zscore, gzscore
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
>>> mu, sigma = 3., 1. # mean and standard deviation
>>> x = rng.lognormal(mu, sigma, size=500)
显示样本的直方图：
>>> fig, ax = plt.subplots()
>>> ax.hist(x, 50)
>>> plt.show()
显示经典 z 分数标准化样本的直方图。分布被重新缩放，但其形状保持不变。
>>> fig, ax = plt.subplots()
>>> ax.hist(zscore(x), 50)
>>> plt.show()
证明几何 z 分数的分布被重新缩放并且准正态：
>>> fig, ax = plt.subplots()
>>> ax.hist(gzscore(x), 50)
>>> plt.show()
scipy.stats.wasserstein_distance
scipy.stats.wasserstein_distance(u_values, v_values, u_weights=None, v_weights=None)
计算两个 1D 分布之间的第一个 Wasserstein 距离。
这个距离也被称为“地球移动距离”，因为它可以看作是将 (u) 转换为 (v) 所需的最小“工作”量，其中“工作”被定义为必须移动的分布权重乘以它必须移动的距离。
新版本 1.0.0 中新增。
参数：
u_values, v_values 数组样式
观察到的值在（经验）分布中。
u_weights, v_weights 数组样式，可选
每个值的权重。如果未指定，每个值将被赋予相同的权重。u_weights（或v_weights）必须与u_values（或v_values）具有相同的长度。如果权重之和与 1 不同，则它仍必须是正的和有限的，以便权重可以归一化为总和为 1。
返回：
distance 浮点数
分布之间计算的距离。
注意事项
分布 (u) 和 (v) 之间的第一个 Wasserstein 距离是：
[l_1 (u, v) = \inf_{\pi \in \Gamma (u, v)} \int_{\mathbb{R} \times \mathbb{R}} |x-y| \mathrm{d} \pi (x, y)]
其中 (\Gamma (u, v)) 是在第一个和第二因子上的边际分布分别为 (u) 和 (v) 的（概率）分布集合。
如果 (U) 和 (V) 是 (u) 和 (v) 的累积分布函数，则此距离也等于：
[l_1(u, v) = \int_{-\infty}^{+\infty} |U-V|]
请参见[2]以证明两个定义的等价性。
输入分布可以是经验的，因此来自样本，其值有效地成为函数的输入，或者它们可以被视为广义函数，此时它们是位于指定值处的 Dirac delta 函数的加权和。
参考文献
[1]
“Wasserstein metric”，en.wikipedia.org/wiki/Wasserstein_metric
[2]
Ramdas, Garcia, Cuturi "On Wasserstein Two Sample Testing and Related Families of Nonparametric Tests" (2015). arXiv:1509.02237.
示例
>>> from scipy.stats import wasserstein_distance
>>> wasserstein_distance([0, 1, 3], [5, 6, 8])
5.0
>>> wasserstein_distance([0, 1], [0, 1], [3, 1], [2, 2])
0.25
>>> wasserstein_distance([3.4, 3.9, 7.5, 7.8], [4.5, 1.4],
... [1.4, 0.9, 3.1, 7.2], [3.2, 3.5])
4.0781331438047861
scipy.stats.energy_distance
scipy.stats.energy_distance(u_values, v_values, u_weights=None, v_weights=None)
计算两个一维分布之间的能量距离。
1.0.0 版新功能。
参数：
u_values, v_valuesarray_like
观察到的（经验）分布中的值。
u_weights, v_weightsarray_like，可选
每个值的权重。如果未指定，则每个值被分配相同的权重。u_weights（v_weights）必须与u_values（v_values）具有相同的长度。如果权重之和不等于 1，则必须仍为正且有限，以便能够将权重归一化为 1。
返回：
distancefloat
计算的分布之间的距离。
注释
两个分布(u)和(v)之间的能量距离，其累积分布函数分别为(U)和(V)，等于：
[D(u, v) = \left( 2\mathbb E|X - Y| - \mathbb E|X - X'| - \mathbb E|Y - Y'| \right)^{1/2}]
其中(X)和(X')（分别(Y)和(Y')）是独立随机变量，其概率分布为(u)（(v)）。
有时，该量的平方被称为“能量距离”（例如在[2]，[4]），但正如[1]和[3]中所指出的那样，仅上述定义符合距离函数（度量）的公理。
如[2]所示，对于一维实值变量，能量距离与 Cramér-von Mises 距离的非分布自由版本相关联：
[D(u, v) = \sqrt{2} l_2(u, v) = \left( 2 \int_{-\infty}^{+\infty} (U-V)² \right)^{1/2}]
注意，普通的 Cramér-von Mises 标准使用距离的无分布版本。详见[2]（第二部分），关于距离两个版本的更多详细信息。
输入分布可以是经验性的，因此来自其值有效作为函数的输入的样本，或者可以视为广义函数，此时它们是位于指定值处的 Dirac δ函数的加权和。
参考文献
[1]
Rizzo, Szekely，“Energy distance.” Wiley Interdisciplinary Reviews: Computational Statistics，8(1):27-38（2015）。
[2] (1,2,3)
Szekely，“E-statistics: The energy of statistical samples.” Bowling Green State University, Department of Mathematics and Statistics, Technical Report 02-16（2002）。
[3]
“Energy distance”，en.wikipedia.org/wiki/Energy_distance
[4]
Bellemare, Danihelka, Dabney, Mohamed, Lakshminarayanan, Hoyer, Munos，“The Cramer Distance as a Solution to Biased Wasserstein Gradients”（2017）。arXiv:1705.10743。
示例
>>> from scipy.stats import energy_distance
>>> energy_distance([0], [2])
2.0000000000000004
>>> energy_distance([0, 8], [0, 8], [3, 1], [2, 2])
1.0000000000000002
>>> energy_distance([0.7, 7.4, 2.4, 6.8], [1.4, 8. ],
... [2.1, 4.2, 7.4, 8. ], [7.6, 8.8])
0.88003340976158217
scipy.stats.rvs_ratio_uniforms
scipy.stats.rvs_ratio_uniforms(pdf, umax, vmin, vmax, size=1, c=0, random_state=None)
使用比例均匀方法从概率密度函数生成随机样本。
自版本 1.12.0 弃用：rvs_ratio_uniforms将在 SciPy 1.15.0 中移除，推荐使用scipy.stats.sampling.RatioUniforms替代。
参数：
pdfcallable
签名为*pdf(x)*的函数，与分布的概率密度函数成比例。
umax浮点数
u-方向边界矩形的上限。
vmin浮点数
v-方向边界矩形的下限。
vmax浮点数
v-方向边界矩形的上限。
size整数或整数元组，可选
定义随机变量的数量（默认为 1）。
c浮点数，可选。
比例均匀方法的偏移参数，请参见注意事项。默认为 0。
random_state{None, int, numpy.random.Generator,
如果seed为 None（或np.random），则使用numpy.random.RandomState单例。如果seed为 int，则使用新的RandomState实例，并使用seed作为种子。如果seed已经是Generator或RandomState实例，则使用该实例。
返回：
rvsndarray
根据概率密度函数定义的随机变量。
注意事项
请参阅scipy.stats.sampling.RatioUniforms获取文档。
scipy.stats.fit
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.fit.html#scipy.stats.fit
scipy.stats.fit(dist, data, bounds=None, *, guess=None, method='mle', optimizer=<function differential_evolution>)
将离散或连续分布拟合到数据
给定分布、数据和分布参数的边界，返回参数的最大似然估计。
参数：
distscipy.stats.rv_continuous 或 scipy.stats.rv_discrete
表示要适应数据的分布对象。
data1D array_like
要进行分布拟合的数据。如果数据包含任何 np.nan、np.inf 或 -np.inf，拟合方法将引发 ValueError。
boundsdict 或 元组序列，可选
如果是字典，则每个键是分布参数的名称，相应的值是该参数的下界和上界组成的元组。如果分布仅在该参数的有限值范围内定义，那么不需要该参数的条目；例如，某些分布的参数必须在 [0, 1] 区间内。位置 (loc) 和尺度 (scale) 参数的边界是可选的；默认情况下，它们分别固定为 0 和 1。
如果是一个序列，第 i 个元素是分布的第 i 个参数的下界和上界组成的元组。在这种情况下，必须提供所有分布形状参数的边界。可选地，位置和尺度的边界可以跟随分布形状参数。
如果要固定某个形状（例如已知的形状），则下界和上界可以相等。如果用户提供的下界或上界超出了分布定义域的边界，分布域的边界将替换用户提供的值。类似地，必须为整数的参数将被限制为用户提供边界内的整数值。
guessdict 或 array_like，可选
如果是字典，则每个键是分布的参数名称，相应的值是参数值的猜测。
如果是一个序列，第 i 个元素是分布的第 i 个参数的猜测值。在这种情况下，必须提供所有分布形状参数的猜测值。
如果未提供 guess，则不会将决策变量的猜测传递给优化器。如果提供了 guess，则会将任何缺失参数的猜测设置为下界和上界的均值。必须为整数的参数的猜测值将四舍五入为整数值，而位于用户提供边界和分布定义域交集之外的猜测值将被剪裁。
method{‘mle’, ‘mse’}
使用method="mle"（默认），通过最小化负对数似然函数来计算拟合。对于超出分布支持的观测值，应用大的有限惩罚（而不是无限的负对数似然）。使用method="mse"，通过最小化负对数产品间距函数来计算拟合。对于超出支持的观测值，应用相同的惩罚。我们遵循[1]的方法，该方法适用于具有重复观测样本。
optimizercallable, optional
optimizer是一个可调用对象，接受以下位置参数。
funcallable
要优化的目标函数。fun接受一个参数x，分布的候选形状参数，并返回给定x、dist和提供的data的目标函数值。optimizer的工作是找到最小化fun的决策变量值。
optimizer还必须接受以下关键字参数。
boundssequence of tuples
决策变量值的边界；每个元素将是包含决策变量下限和上限的元组。
如果提供了guess，optimizer还必须接受以下关键字参数。
x0array_like
每个决策变量的猜测。
如果分布有任何必须是整数的形状参数，或者如果分布是离散的且位置参数不固定，optimizer还必须接受以下关键字参数。
integralityarray_like of bools
对于每个决策变量，如果决策变量必须被限制为整数值，则为 True，如果决策变量是连续的，则为 False。
optimizer必须返回一个对象，例如scipy.optimize.OptimizeResult的实例，其中将决策变量的最优值保存在属性x中。如果提供了fun、status或message属性，它们将包含在由fit返回的结果对象中。
返回：
resultFitResult
具有以下字段的对象。
paramsnamedtuple
包含分布形状参数、位置和（如果适用）尺度的最大似然估计的命名元组。
successbool 或 None
优化器是否考虑优化是否成功终止。
messagestr 或 None
优化器提供的任何状态消息。
对象具有以下方法：
nllf(params=None, data=None)
默认情况下，负对数似然函数在给定数据的拟合params处。接受包含分布的替代形状、位置和尺度以及替代数据数组的元组。
plot(ax=None)
将拟合分布的 PDF/PMF 叠加在数据的归一化直方图上。
另请参见
rv_continuous, rv_discrete
注意事项
当用户提供包含最大似然估计的紧密界限时，优化更有可能收敛到最大似然估计。例如，当将二项分布拟合到数据时，每个样本背后的实验数可能已知，在这种情况下，相应的形状参数n可以固定。
参考文献
[1]
邵勇照和 Marjorie G. Hahn. “最大间隔乘积方法：具有强一致性的统一表达。” 伊利诺伊数学期刊 43.3 (1999): 489-499。
例子
假设我们希望将分布拟合到以下数据。
>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> dist = stats.nbinom
>>> shapes = (5, 0.5)
>>> data = dist.rvs(*shapes, size=1000, random_state=rng)
假设我们不知道数据是如何生成的，但我们怀疑它遵循负二项分布，参数为n和p。 (参见scipy.stats.nbinom.) 我们相信参数n小于 30，且参数p必须在区间[0, 1]内。我们将这些信息记录在变量bounds中，并将其传递给fit。
>>> bounds = [(0, 30), (0, 1)]
>>> res = stats.fit(dist, data, bounds)
fit在用户指定的bounds范围内搜索最佳与数据匹配的值（以最大似然估计的意义）。在这种情况下，它找到了与实际生成数据相似的形状值。
>>> res.params
FitParams(n=5.0, p=0.5028157644634368, loc=0.0) # may vary
我们可以通过在数据的归一化直方图上叠加分布的概率质量函数（形状适合数据）来可视化结果。
>>> import matplotlib.pyplot as plt # matplotlib must be installed to plot
>>> res.plot()
>>> plt.show()
注意，n的估计值恰好是整数；这是因为nbinom PMF 的定义域仅包含整数n，而nbinom对象“知道”这一点。nbinom还知道，形状p必须是介于 0 和 1 之间的值。在这种情况下 - 当分布的域对于参数是有限的时候 - 我们不需要为参数指定边界。
>>> bounds = {'n': (0, 30)} # omit parameter p using a `dict`
>>> res2 = stats.fit(dist, data, bounds)
>>> res2.params
FitParams(n=5.0, p=0.5016492009232932, loc=0.0) # may vary
如果我们希望强制分布在n固定为 6 的情况下进行拟合，我们可以将n的上下界都设为 6。然而，请注意，在这种情况下，优化的目标函数值通常会更差（更高）。
>>> bounds = {'n': (6, 6)} # fix parameter `n`
>>> res3 = stats.fit(dist, data, bounds)
>>> res3.params
FitParams(n=6.0, p=0.5486556076755706, loc=0.0) # may vary
>>> res3.nllf() > res.nllf()
True # may vary
注意，前述示例的数值结果是典型的，但可能会有所不同，因为fit使用的默认优化器scipy.optimize.differential_evolution是随机的。然而，我们可以通过定制优化器的设置来确保可复现性 - 或者完全使用不同的优化器 - 使用optimizer参数。
>>> from scipy.optimize import differential_evolution
>>> rng = np.random.default_rng()
>>> def optimizer(fun, bounds, *, integrality):
... return differential_evolution(fun, bounds, strategy='best2bin',
... seed=rng, integrality=integrality)
>>> bounds = [(0, 30), (0, 1)]
>>> res4 = stats.fit(dist, data, bounds, optimizer=optimizer)
>>> res4.params
FitParams(n=5.0, p=0.5015183149259951, loc=0.0)
scipy.stats.ecdf
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.ecdf.html#scipy.stats.ecdf
scipy.stats.ecdf(sample)
样本的经验累积分布函数。
经验累积分布函数（ECDF）是样本底层分布的 CDF 的阶梯函数估计。此函数返回表示经验分布函数及其补集经验生存函数的对象。
参数：
sample1D 数组或scipy.stats.CensoredData
除了数组，支持scipy.stats.CensoredData的实例，包括未审查和右审查的观察结果。目前，其他scipy.stats.CensoredData的实例将导致NotImplementedError。
返回：
resECDFResult
具有以下属性的对象。
cdfEmpiricalDistributionFunction
表示样本的经验累积分布函数的对象。
sfEmpiricalDistributionFunction
表示经验生存函数的对象。
cdf和sf属性本身具有以下属性。
quantilesndarray
定义经验 CDF/SF 的样本中的唯一值。
probabilitiesndarray
指数的分位数对应的概率点估计。
和以下方法：
evaluate(x)：
在参数处评估 CDF/SF。
plot(ax)：
在提供的坐标轴上绘制 CDF/SF。
confidence_interval(confidence_level=0.95)：
计算在分位数值周围 CDF/SF 的置信区间。
注意事项
当样本的每个观测是精确测量时，ECDF 在每个观测点[1]处按1/len(sample)递增。
当观测值为下限、上限或上下限时，数据被称为“审查”，sample可以作为scipy.stats.CensoredData的实例提供。
对于右审查数据，ECDF 由 Kaplan-Meier 估计器给出[2]；目前不支持其他形式的审查。
置信区间根据格林伍德公式或更近期的“指数格林伍德”公式计算，如[4]所述。
参考文献
[1] (1,2,3)
Conover, William Jay. 实用非参数统计. Vol. 350. John Wiley & Sons, 1999.
[2]
Kaplan, Edward L., and Paul Meier. “非参数估计来自不完整观测。” 美国统计协会杂志 53.282 (1958): 457-481.
[3]
Goel, Manish Kumar, Pardeep Khanna, and Jugal Kishore. “理解生存分析：Kaplan-Meier 估计。” 国际阿育吠陀研究杂志 1.4 (2010): 274.
[4]
Sawyer, Stanley. “生存分析中的格林伍德和指数格林伍德置信区间。” www.math.wustl.edu/~sawyer/handouts/greenwood.pdf
示例
非截尾数据
如在示例[1]第 79 页中，从一所高中的男生中随机选择了五个男生。他们的一英里跑步时间记录如下。
>>> sample = [6.23, 5.58, 7.06, 6.42, 5.20] # one-mile run times (minutes)
经验分布函数，用来近似样本男生一英里跑步时间的总体分布函数，计算如下。
>>> from scipy import stats
>>> res = stats.ecdf(sample)
>>> res.cdf.quantiles
array([5.2 , 5.58, 6.23, 6.42, 7.06])
>>> res.cdf.probabilities
array([0.2, 0.4, 0.6, 0.8, 1\. ])
要将结果绘制为阶跃函数：
>>> import matplotlib.pyplot as plt
>>> ax = plt.subplot()
>>> res.cdf.plot(ax)
>>> ax.set_xlabel('One-Mile Run Time (minutes)')
>>> ax.set_ylabel('Empirical CDF')
>>> plt.show()
右截尾数据
如在示例[1]第 91 页中，对十个汽车传动带的使用寿命进行了测试。五次测试因测试中的传动带断裂而结束，但其余测试因其他原因结束（例如，研究资金耗尽，但传动带仍然功能良好）。记录了传动带的行驶里程如下。
>>> broken = [77, 47, 81, 56, 80] # in thousands of miles driven
>>> unbroken = [62, 60, 43, 71, 37]
在测试结束时仍然功能良好的传动带的精确寿命时间是未知的，但已知它们超过了记录在unbroken中的值。因此，这些观察被称为“右截尾”，并且使用scipy.stats.CensoredData来表示数据。
>>> sample = stats.CensoredData(uncensored=broken, right=unbroken)
经验生存函数计算如下。
>>> res = stats.ecdf(sample)
>>> res.sf.quantiles
array([37., 43., 47., 56., 60., 62., 71., 77., 80., 81.])
>>> res.sf.probabilities
array([1\. , 1\. , 0.875, 0.75 , 0.75 , 0.75 , 0.75 , 0.5 , 0.25 , 0\. ])
要将结果绘制为阶跃函数：
>>> ax = plt.subplot()
>>> res.cdf.plot(ax)
>>> ax.set_xlabel('Fanbelt Survival Time (thousands of miles)')
>>> ax.set_ylabel('Empirical SF')
>>> plt.show()
scipy.stats.logrank
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.logrank.html#scipy.stats.logrank
scipy.stats.logrank(x, y, alternative='two-sided')
通过 logrank 测试比较两个样本的生存分布。
参数：
x, yarray_like or CensoredData
根据其经验生存函数比较样本。
alternative{‘two-sided’, ‘less’, ‘greater’}, optional
定义备择假设。
零假设是两组，如X和Y的生存分布相同。
可用以下备择假设[4]（默认为‘two-sided’）：
-
‘two-sided’: 两组的生存分布不相同。
-
‘less’: 生存组X更受青睐：在某些时间点上，组X的失效率函数小于组Y的失效率函数。
-
‘greater’: survival of group Y is favored: the group X failure rate function is greater than the group Y failure rate function at some times.
返回：
resLogRankResult
一个包含属性的对象：
statisticfloat ndarray
所计算的统计量（如下定义）。其大小是多数其他 logrank 测试实现返回的平方根大小。
pvaluefloat ndarray
测试的计算 p 值。
参见
scipy.stats.ecdf
注意事项
logrank 测试[1]比较观察到的事件数与零假设下预期事件数之间的差异，即两个样本是否从相同分布中抽取。统计量为
[Z_i = \frac{\sum_{j=1}^J(O_{i,j}-E_{i,j})}{\sqrt{\sum_{j=1}^J V_{i,j}}} \rightarrow \mathcal{N}(0,1)]
where
[E_{i,j} = O_j \frac{N_{i,j}}{N_j}, \qquad V_{i,j} = E_{i,j} \left(\frac{N_j-O_j}{N_j}\right) \left(\frac{N_j-N_{i,j}}{N_j-1}\right),]
(i) denotes the group (i.e. it may assume values (x) or (y), or it may be omitted to refer to the combined sample) (j) denotes the time (at which an event occurred), (N) is the number of subjects at risk just before an event occurred, and (O) is the observed number of events at that time.
由logrank返回的statistic (Z_x)是许多其他实现返回的统计量的（带符号的）平方根。在零假设下，(Z_x**2)渐近地按自由度为一的卡方分布分布。因此，(Z_x)渐近地按标准正态分布分布。使用(Z_x)的优势在于保留了符号信息（即观察到的事件数是否倾向于少于或大于零假设下预期的数量），从而允许scipy.stats.logrank提供单侧备择假设。
参考文献
[1]
Mantel N. “评估生存数据及其相关的两个新秩次统计量。”《癌症化疗报告》，50(3):163-170，PMID: 5910392，1966 年
[2]
Bland, Altman, “对数秩检验”，BMJ，328:1073，DOI:10.1136/bmj.328.7447.1073，2004 年
[3]
“对数秩检验”，维基百科，zh.wikipedia.org/wiki/对数秩检验
[4]
Brown, Mark. “关于对数秩检验方差选择的问题。”《生物统计学》，71.1 (1984): 65-74.
[5]
Klein, John P., 和 Melvin L. Moeschberger.《生存分析：截尾和删节数据的技术》。卷 1230. 纽约：Springer，2003 年。
示例
参考文献[2] 比较了两种不同类型复发性恶性胶质瘤患者的生存时间。下面的样本记录了每位患者参与研究的时间（以周为单位）。由于数据是右截尾的：未截尾的观察对应于观察到的死亡，而截尾的观察对应于患者因其他原因离开研究，因此使用了scipy.stats.CensoredData 类。
>>> from scipy import stats
>>> x = stats.CensoredData(
... uncensored=[6, 13, 21, 30, 37, 38, 49, 50,
... 63, 79, 86, 98, 202, 219],
... right=[31, 47, 80, 82, 82, 149]
... )
>>> y = stats.CensoredData(
... uncensored=[10, 10, 12, 13, 14, 15, 16, 17, 18, 20, 24, 24,
... 25, 28,30, 33, 35, 37, 40, 40, 46, 48, 76, 81,
... 82, 91, 112, 181],
... right=[34, 40, 70]
... )
我们可以计算和可视化两组的经验生存函数如下。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> ax = plt.subplot()
>>> ecdf_x = stats.ecdf(x)
>>> ecdf_x.sf.plot(ax, label='Astrocytoma')
>>> ecdf_y = stats.ecdf(y)
>>> ecdf_x.sf.plot(ax, label='Glioblastoma')
>>> ax.set_xlabel('Time to death (weeks)')
>>> ax.set_ylabel('Empirical SF')
>>> plt.legend()
>>> plt.show()
经验生存函数的视觉检查表明，两组的生存时间倾向于不同。为了正式评估这种差异是否在 1%水平上显著，我们使用了对数秩检验。
>>> res = stats.logrank(x=x, y=y)
>>> res.statistic
-2.73799...
>>> res.pvalue
0.00618...
p 值小于 1%，因此我们可以认为数据证据不支持零假设，支持备择假设，即两个生存函数之间存在差异。
scipy.stats.directional_stats
scipy.stats.directional_stats(samples, *, axis=0, normalize=True)
计算方向数据的样本统计量。
计算样本向量的方向平均值（也称为平均方向向量）和平均结果长度。
方向平均值是向量数据的“首选方向”的度量。它类似于样本均值，但在数据长度无关紧要时使用（例如单位向量）。
平均结果长度是一个介于 0 和 1 之间的值，用于量化方向数据的分散程度：平均结果长度越小，分散程度越大。关于涉及平均结果长度的方向方差的多个定义可见[1]和[2]。
参数：
samplesarray_like
输入数组。必须至少是二维的，且输入的最后一个轴必须与向量空间的维数对应。当输入恰好是二维时，这意味着数据的每一行都是一个向量观测值。
axis整数，默认为 0
计算方向平均值的轴。
normalize: 布尔值，默认为 True
如果为 True，则将输入标准化，以确保每个观测值都是单位向量。如果观测值已经是单位向量，则考虑将其设置为 False，以避免不必要的计算。
返回：
resDirectionalStats
一个包含属性的对象：
mean_directionndarray
方向平均值。
mean_resultant_lengthndarray
平均结果长度 [1]。
另请参阅
circmean
循环均值；即 2D 角度的方向均值。
circvar
循环方差；即 2D 角度的方向方差。
注意事项
此处使用了来自[1]的方向平均值定义。假设观测值是单位向量，则计算如下。
mean = samples.mean(axis=0)
mean_resultant_length = np.linalg.norm(mean)
mean_direction = mean / mean_resultant_length
此定义适用于方向数据（即每个观测的大小无关紧要的向量数据），但不适用于轴向数据（即每个观测的大小和符号都无关紧要的向量数据）。
已经提出了几个涉及平均结果长度 R 的方向方差的定义，包括 1 - R [1]、 1 - R**2 [2] 和 2 * (1 - R) [2]。与选择其中一个不同，此函数返回 R 作为属性 mean_resultant_length，以便用户可以计算其首选的分散度量。
参考
[1] (1,2,3,4)
Mardia, Jupp. (2000). Directional Statistics (p. 163). Wiley.
[2] (1,2,3,4)
en.wikipedia.org/wiki/Directional_statistics
示例
>>> import numpy as np
>>> from scipy.stats import directional_stats
>>> data = np.array([[3, 4], # first observation, 2D vector space
... [6, -8]]) # second observation
>>> dirstats = directional_stats(data)
>>> dirstats.mean_direction
array([1., 0.])
相比之下，向量的常规样本均值会受每个观测值的大小的影响。此外，结果不会是一个单位向量。
>>> data.mean(axis=0)
array([4.5, -2.])
directional_stats的一个典型用例是在球面上寻找一组观测值的有意义中心，例如地理位置。
>>> data = np.array([[0.8660254, 0.5, 0.],
... [0.8660254, -0.5, 0.]])
>>> dirstats = directional_stats(data)
>>> dirstats.mean_direction
array([1., 0., 0.])
另一方面，常规样本均值的结果不位于球面表面上。
>>> data.mean(axis=0)
array([0.8660254, 0., 0.])
该函数还返回平均结果长度，可用于计算方向方差。例如，使用定义 Var(z) = 1 - R 来自于[2]，其中 R 是平均结果长度，我们可以计算上述示例中向量的方向方差为：
>>> 1 - dirstats.mean_resultant_length
0.13397459716167093
scipy.stats.circmean
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.circmean.html#scipy.stats.circmean
scipy.stats.circmean(samples, high=6.283185307179586, low=0, axis=None, nan_policy='propagate', *, keepdims=False)
计算范围内样本的圆均值。
参数：
samples类似数组
输入数组。
highfloat 或 int，可选
样本范围的高边界。默认为 2*pi。
lowfloat 或 int，可选
样本范围的低边界。默认为 0。
axisint 或 None，默认为 None
如果是 int，则为计算统计量的输入轴。输入的每个轴切片（例如行）的统计量将出现在输出的相应元素中。如果为 None，则在计算统计量之前将对输入进行展平。
nan_policy{'propagate', 'omit', 'raise'}
定义如何处理输入的 NaN。
-
propagate: 如果计算统计量的轴切片（例如行）中存在 NaN，则输出的相应条目将为 NaN。 -
omit: 在执行计算时将忽略 NaN。如果在计算统计量的轴切片中剩余的数据不足，则输出的相应条目将为 NaN。 -
raise: 如果存在 NaN，则会引发ValueError。
keepdimsbool，默认为 False
如果设置为 True，则减少的轴将作为尺寸为一的维度保留在结果中。选择此选项可确保结果正确地广播到输入数组。
返回：
circmeanfloat
圆均值。
另请参阅
circstd
圆标准差。
circvar
圆方差。
注意事项
从 SciPy 1.9 开始，np.matrix 输入（不建议用于新代码）在执行计算之前将转换为 np.ndarray。在这种情况下，输出将是适当形状的标量或 np.ndarray，而不是二维的 np.matrix。同样，虽然忽略掩码数组的屏蔽元素，输出将是标量或 np.ndarray，而不是具有 mask=False 的掩码数组。
示例
为简单起见，所有角度都以度数打印出来。
>>> import numpy as np
>>> from scipy.stats import circmean
>>> import matplotlib.pyplot as plt
>>> angles = np.deg2rad(np.array([20, 30, 330]))
>>> circmean = circmean(angles)
>>> np.rad2deg(circmean)
7.294976657784009
>>> mean = angles.mean()
>>> np.rad2deg(mean)
126.66666666666666
绘制并比较圆均值与算术平均值。
>>> plt.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
... np.sin(np.linspace(0, 2*np.pi, 500)),
... c='k')
>>> plt.scatter(np.cos(angles), np.sin(angles), c='k')
>>> plt.scatter(np.cos(circmean), np.sin(circmean), c='b',
... label='circmean')
>>> plt.scatter(np.cos(mean), np.sin(mean), c='r', label='mean')
>>> plt.legend()
>>> plt.axis('equal')
>>> plt.show()
scipy.stats.circvar
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.stats.circvar.html#scipy.stats.circvar
scipy.stats.circvar(samples, high=6.283185307179586, low=0, axis=None, nan_policy='propagate', *, keepdims=False)
计算假定在范围内的样本的圆形方差。
参数：
samplesarray_like
输入数组。
highfloat 或 int，可选
样本范围的高边界。默认为 2*pi。
lowfloat 或 int，可选
样本范围的低边界。默认为 0。
axisint 或 None，默认为 None
如果是 int，则是计算统计量的输入的轴。输入的每个轴切片（例如行）的统计量将出现在输出的相应元素中。如果是 None，则在计算统计量之前会将输入展平。
nan_policy{‘propagate’, ‘omit’, ‘raise’}
定义如何处理输入的 NaN 值。
-
propagate: 如果在计算统计量的轴切片（例如行）中存在 NaN，则输出的相应条目将是 NaN。 -
omit: 在执行计算时将省略 NaN。如果在计算统计量的轴切片中剩余的数据不足，则输出的相应条目将为 NaN。 -
raise: 如果存在 NaN，则会引发ValueError。
keepdimsbool，默认为 False
如果设置为 True，则被缩减的轴会以尺寸为一的维度保留在结果中。使用此选项，结果将正确地与输入数组进行广播。
返回：
circvarfloat
圆形方差。
参见
circmean
圆形平均值。
circstd
圆形标准差。
注意
这里使用的圆形方差的定义是 1-R，其中 R 是平均结果向量。返回的值在范围 [0, 1] 内，0 表示无方差，1 表示大方差。在小角度的极限情况下，该值类似于“线性”方差的一半。
从 SciPy 1.9 开始，np.matrix 输入（不建议新代码使用）在执行计算之前会转换为 np.ndarray。在这种情况下，输出将是一个适当形状的标量或 np.ndarray，而不是一个 2D 的 np.matrix。同样地，忽略掩码数组的掩码元素后，输出将是一个标量或 np.ndarray，而不是带有 mask=False 的掩码数组。
参考文献
[1]
Fisher, N.I. Statistical analysis of circular data. Cambridge University Press, 1993.
示例
>>> import numpy as np
>>> from scipy.stats import circvar
>>> import matplotlib.pyplot as plt
>>> samples_1 = np.array([0.072, -0.158, 0.077, 0.108, 0.286,
... 0.133, -0.473, -0.001, -0.348, 0.131])
>>> samples_2 = np.array([0.111, -0.879, 0.078, 0.733, 0.421,
... 0.104, -0.136, -0.867, 0.012, 0.105])
>>> circvar_1 = circvar(samples_1)
>>> circvar_2 = circvar(samples_2)
绘制样本。
>>> fig, (left, right) = plt.subplots(ncols=2)
>>> for image in (left, right):
... image.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
... np.sin(np.linspace(0, 2*np.pi, 500)),
... c='k')
... image.axis('equal')
... image.axis('off')
>>> left.scatter(np.cos(samples_1), np.sin(samples_1), c='k', s=15)
>>> left.set_title(f"circular variance: {np.round(circvar_1, 2)!r}")
>>> right.scatter(np.cos(samples_2), np.sin(samples_2), c='k', s=15)
>>> right.set_title(f"circular variance: {np.round(circvar_2, 2)!r}")
>>> plt.show()
