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, rng=None)[源代码]#

执行拟合优度检验,比较数据与分布族。

给定一个分布族和数据,执行零假设检验,即数据是从该族中的分布中抽取的。可以指定分布的任何已知参数。分布的剩余参数将拟合到数据,并相应地计算检验的 p 值。有几种统计量可用于比较分布与数据。

参数:
distscipy.stats.rv_continuous

表示零假设下的分布族的对象。

data一维类数组

要测试的有限的、非截尾的数据。

known_paramsdict,可选

一个字典,包含已知分布参数的名称-值对。蒙特卡罗样本是从具有这些参数值的零假设分布中随机抽取的。在评估观察到的 data 和每个蒙特卡罗样本的统计量之前,仅将零假设分布族的剩余未知参数拟合到样本;已知参数保持固定。如果分布族的所有参数都是已知的,则省略将分布族拟合到每个样本的步骤。

fit_paramsdict,可选

一个字典,包含已拟合到数据的分布参数的名称-值对,例如使用 scipy.stats.fitdistfit 方法。蒙特卡罗样本是从具有这些指定参数值的零假设分布中抽取的。但是,无论样本是观察到的 data 还是蒙特卡罗样本,这些参数和零假设分布族的所有其他未知参数始终拟合到样本,然后才评估统计量。

guessed_paramsdict,可选

一个字典,包含已猜测的分布参数的名称-值对。这些参数始终被视为自由参数,并且既拟合到提供的 data,也拟合到从零假设分布中抽取的蒙特卡罗样本。这些 guessed_params 的目的是用作数值拟合过程的初始值。

statistic{“ad”, “ks”, “cvm”, “filliben”} 或可调用对象,可选

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

n_mc_samplesint,默认值:9999

从零假设分布中抽取的蒙特卡罗样本的数量,以形成统计量的零分布。每个样本的大小与给定的 data 相同。

rng{None, int, numpy.random.Generator},可选

如果通过关键字传递 rng,则将 numpy.random.Generator 之外的类型传递给 numpy.random.default_rng 以实例化 Generator。如果 rng 已经是 Generator 实例,则使用提供的实例。指定 rng 以获得可重复的函数行为。

如果此参数是通过位置传递的,或者 random_state 是通过关键字传递的,则应用参数 random_state 的旧行为

  • 如果 random_state 为 None(或 numpy.random),则使用 numpy.random.RandomState 单例。

  • 如果 random_state 是一个整数,则使用一个新的 RandomState 实例,并使用 random_state 作为种子。

  • 如果 random_state 已经是 GeneratorRandomState 实例,则使用该实例。

在版本 1.15.0 中更改:作为 SPEC-007 从使用 numpy.random.RandomState 转换为 numpy.random.Generator 的一部分,此关键字已从 random_state 更改为 rng。在过渡期内,这两个关键字将继续工作,但一次只能指定一个。在过渡期之后,使用 random_state 关键字的函数调用将发出警告。random_staterng 的行为如上所述,但在新代码中应仅使用 rng 关键字。

返回:
resGoodnessOfFitResult

具有以下属性的对象。

fit_resultFitResult

表示提供的 dist 拟合到 data 的对象。此对象包括完全定义零假设分布(即从中抽取蒙特卡罗样本的分布)的分布族参数的值。

statisticfloat

比较提供的 data 与零假设分布的统计量的值。

pvaluefloat

零分布中统计量值至少与提供的 data 的统计量值一样极端的元素的比例。

null_distributionndarray

从零假设分布中抽取的每个蒙特卡罗样本的统计量值。

注意

这是一个广义的蒙特卡罗拟合优度过程,其特殊情况与各种安德森-达林检验、Lilliefors 检验等相对应。该检验在 [2][3][4] 中描述为参数自举检验。这是一个蒙特卡罗检验,其中指定从中抽取样本的分布的参数已从数据中估计出来。我们自始至终使用“蒙特卡罗”而不是“参数自举”来描述检验,以避免与更熟悉的非参数自举混淆,并描述如下如何执行检验。

传统的拟合优度检验

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

算法概述

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

  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\) 是蒙特卡洛零分布中大于或等于为数据计算的统计量值的统计量值的数量,\(m\) 是蒙特卡洛零分布中元素的数量(n_mc_samples)。分子和分母加 \(1\) 可以被认为是将与数据对应的统计量值包含在零分布中,但在[5]中给出了更正式的解释。

局限性

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

反模式

因此,可能很想将预先拟合到数据(由用户)的分布参数视为已知参数,因为指定分布的所有参数排除了将分布拟合到每个蒙特卡洛样本的需要。(这基本上是执行原始 Kilmogorov-Smirnov 测试的方式。)尽管此类测试可以提供反对零假设的证据,但该测试是保守的,因为小的 p 值往往会(大大)高估犯 I 类错误(即,尽管零假设为真,但拒绝零假设)的概率,并且该测试的功效很低(即,即使零假设为假,它也不太可能拒绝零假设)。这是因为蒙特卡洛样本不太可能像数据那样与零假设分布一致。这往往会增加零分布中记录的统计量值,从而使其中较大一部分超过数据的统计量值,从而夸大 p 值。

参考文献

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

M. A. Stephens (1974)。“拟合优度和一些比较的 EDF 统计量。”美国统计协会杂志,第 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)。“半参数模型中用于拟合优度检验的参数自举的有效性。”Annales de l’IHP Probabilités et statistiques。第 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)。“关于均值和方差未知的正态性的 Kolmogorov-Smirnov 检验。”美国统计协会杂志 62.318: 399-402。

[7]

Filliben,James J。“正态性的概率图相关系数检验。”Technometrics 17.1 (1975): 111-117。

示例

一个著名的检验数据是否从给定分布中抽样的零假设的检验是 Kolmogorov-Smirnov (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', rng=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 进行比较。这往往会降低观测数据的统计量值,但是在计算其他样本(例如我们随机抽取以形成蒙特卡洛零分布的样本)的统计量时,这是“不公平的”。很容易纠正这种情况:每当我们计算样本的 KS 统计量时,我们都使用拟合到该样本的正态分布的 CDF。在这种情况下,零分布尚未被精确计算,并且通常使用如上所述的蒙特卡洛方法来近似。这就是 goodness_of_fit 的优势所在。

>>> res = stats.goodness_of_fit(stats.norm, x, statistic='ks',
...                             rng=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',
...                             rng=rng)
>>> res.statistic, res.pvalue
(1.2139573337497467, 0.0034)

使用 goodness_of_fit 的另一个优势是,它不局限于特定的分布集或参数是已知还是必须从数据中估计的条件。相反,对于任何具有足够快和可靠的 fit 方法的分布,goodness_of_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}, rng=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)