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{"ad","ks","cvm","filliben"} 或 callable,可选

在拟合分布族的未知参数到数据后用于比较数据和分布的统计信息。Anderson-Darling(“ad”) [1]、Kolmogorov-Smirnov(“ks”) [1]、Cramer-von Mises(“cvm”) [1] 和 Filliben(“filliben”) [7] 统计信息可供使用。也可以提供带有签名 (dist, data, axis) 的可调用项以计算统计信息。此处 dist 是冻结分布对象(可能带有数组参数),data 是蒙特卡罗样本的数组(具有兼容的形状),axis 是必须在其中计算统计信息的 data 的轴。

n_mc_samplesint,默认值:9999

从零假设分布绘制的蒙特卡罗样本数,以形成统计信息的零分布。每个样本的大小与给定的 data 相同。

random_state{None,int,numpy.random.Generator

生成蒙特卡罗样本时使用的伪随机数生成器状态。

如果 random_stateNone(默认值),则使用 numpy.random.RandomState 单例。如果 random_state 是 int,则使用新的 RandomState 实例,用 random_state 播种。如果 random_state 已经是 GeneratorRandomState 实例,则将使用提供的实例。

返回:
resGoodnessOfFitResult

具有以下属性的对象。

fit_resultFitResult

一个表示所提供的 distdata 拟合的对象。此对象包括完全定义原假设分布的分布族参数的值,即从中抽取蒙特卡罗样本的分布。

statisticfloat

用于比较所提供的 data 与原假设分布的统计量值。

pvaluefloat

statistic 值至少与所提供的 data 的 statistic 值一样极端的原假设分布中元素的比例。

null_distributionndarray

从原假设分布中抽取的每个蒙特卡罗样本针对 statistic 的值。

注释

这是一种广义的蒙特卡罗拟合优度程序,其特例对应于各种 Anderson-Darling 检验、Lilliefors 检验等。该检验在 [2][3][4] 中被描述为一种参数引导检验。这是一种蒙特卡罗检验,其中从数据中估算出了用于指定从中抽取样本的分布的参数。我们一直在使用“蒙特卡罗”而非“参数引导”描述此检验,以避免与更常见的非参数引导混淆,并在下面描述如何执行检验。

传统的拟合优度检验

传统上,使用蒙特卡罗方法预先计算出与一组固定显著性水平相对应的临界值。用户仅针对其观测到的 data 计算检验统计量的值并将其与制表临界值进行比较来执行检验。这种做法并不十分灵活,因为无法为所有分布和已知与未知参数值的组合提供表格。此外,当临界值从有限的制表数据中内插以对应于用户的样本量和拟合参数值时,结果可能会不准确。为了克服这些缺点,此函数允许用户针对其特定数据执行蒙特卡罗试验。

算法概述

简而言之,此例程执行以下步骤

  1. 对给定的 data 拟合未知参数,从而形成“原假设”分布,并计算此对数据和分布的 statistic。

  2. 从该原假设分布中抽取随机样本。

  3. 对每个随机样本拟合未知参数。

  4. 计算每个样本与其已拟合到该样本的分布之间的 statistic。

  5. 将对应于 (1) 中数据的统计值与对应于 (4) 中的随机样本的统计值进行比较。p 值是统计值大于或等于观测数据的统计值的样本比例。

更详细地说,步骤如下。

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

接下来,从零假设分布中抽取了许多(具体来说n_mc_samples)新样本,每个样本包含与data相同数量的观测值。每个重新抽样都要拟合分配族dist的所有未知参数,并计算每个样本与其对应的拟合分布之间的statistic。这些统计值形成了蒙特卡罗零分布(不要与上文中的“零假设分布”混淆)。

此测试的 p 值是蒙特卡罗零分布中至少与所提供data的统计值一样极端的统计值的比例。更确切地说,p 值由以下公式得出:

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

其中\(b\)是蒙特卡罗零分布中大于或等于为data计算的统计值的统计值的个数,\(m\)是蒙特卡罗零分布中的元素个数(n_mc_samples)。分子分母加上\(1\)可以认为在零分布中包括了对应于data的统计值,但[5]中给出了更正式的解释。

局限性

对于某些分布族,此测试可能非常慢,因为必须为每个蒙特卡罗样本拟合分布族的未知参数,并且对于 SciPy 中的大多数分布,分布拟合是通过数值优化执行的。

反模式

由于这个原因,将分布参数预先拟合到用户(数据)就可能诱使将它们视为known_params,因为指定分布的所有参数就排除了将分布拟合到每个蒙特卡罗样本中的需要。(这基本上是原始基尔莫格洛夫-斯米尔诺夫检验执行的方式。)尽管这样的检验可以提供反对零假设的证据,但从这个意义上讲,该检验是保守的,即较小的 p 值将倾向于大大高估 I 类错误(即使零假设为真也会拒绝)的可能性,并且检验的功效很低(即使零假设为假,也可能驳回零假设)。这是因为蒙特卡罗样不太可能与零假设分布以及数据一样。这往往会增加在零分布中记录的统计量的值,以至于其中有更多的值会高于数据的统计量值,从而膨胀 p 值。

参考资料

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

M. A. 史蒂文斯(1974)。“拟合优度和一些比较的 EDF 统计数据。”《美国统计协会杂志》第 69 卷,第 730-737 页。

[2]

W. 斯图特、W.G. 曼特加和 M.P. 昆迪米尔(1993)。“基于引导程序的拟合优度检验。”《梅特里卡》40.1:243-256。

[3]

C. 热内斯特和 B. 雷米拉德。(2008)。“半参数模型中拟合优度检验的参数引导程序的有效性。”《法国数学与信息学院纪要概率和统计》第 44 卷,第 6 期。

[4]

I. 科贾迪诺维奇和 J. 颜(2012)。“基于加权引导程序的拟合优度检验:参数引导程序的快速大样本替代方法。”《加拿大统计杂志》40.3:480-500。

[5]

B. 菲普森和 G.K. 史密斯(2010)。“排列 P 值永远不应为零:在随机抽取排列时计算精确 P 值。”《统计在遗传学和分子生物学中的应用》9.1。

[6]

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

[7]

菲利本,詹姆斯 J.“正态分布概率图相关系数检验。”《泰克诺梅特里卡》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,
             statistic_location=0.7751845155861765,
             statistic_sign=-1)

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 检验中,此项统计数据总是将数据与拟合给观测数据的正态分布的 CDF 进行比较。这往往会降低观测数据统计数据的值,但当为其他样本(例如我们将随机提取以形成 Monte Carlo 零分布的样本)计算统计数据时,这种比较“不公平”。对此进行修正很容易:每当我们计算样本的 KS 统计数据时,我们便使用拟合给该样本的正态分布的 CDF。在这种情况下,零分布尚未精确计算,而是使用如上所述的 Monte Carlo 方法进行近似估算。这就是 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 统计数据的最初优势在于能够理论计算零分布,但现在我们能够通过计算近似零分布的方式来使用敏感性更高的统计数据,这将导致更高的检验功效。安德森-达令统计数据 [1] 往往更敏感,并且已经使用 Monte Carlo 方法为各种显着性水平和样本大小编制了该统计数据临界值表。

>>> 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 统计数据针对具有已知位置和未知尺度的瑞利分布执行拟合优度检验。

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

如果分布不完全适合观察到的数据,测试可能无法控制 1 型错误率,即在零假设为真时拒绝零假设的可能性。

我们还应该寻找零分布中的极端离群值,这些极端离群值可能是由不可靠的拟合引起的。这些并不一定使结果无效,但它们往往会降低测试的能力。

>>> _, 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

这个图看起来令人安心。

如果class="docutils literal notranslate">fit方法工作可靠,并且检验统计量的分布对拟合参数的值不是特别敏感,那么goodness_of_fit提供的 p 值有望成为一个良好的近似值。

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