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)[source]#
执行拟合优度检验,比较数据与分布族。
给定一个分布族和数据,执行零假设的检验,即数据是从该族中的一个分布中抽取的。可以指定分布的任何已知参数。分布的剩余参数将拟合到数据,并相应地计算检验的 p 值。有几种统计量可用于比较分布与数据。
- 参数:
- dist
scipy.stats.rv_continuous
表示零假设下的分布族的对象。
- data一维类数组
要测试的有限、非删失数据。
- known_paramsdict,可选
包含已知分布参数的名称-值对的字典。蒙特卡罗样本是从零假设分布中随机抽取的,参数具有这些值。在评估观测 data 和每个蒙特卡罗样本的统计量之前,仅将零假设分布族的剩余未知参数拟合到样本;已知参数保持固定。如果分布族的所有参数都是已知的,那么将分布族拟合到每个样本的步骤将被省略。
- fit_paramsdict,可选
包含已拟合到数据的分布参数的名称-值对的字典,例如使用
scipy.stats.fit
或 dist 的fit
方法。蒙特卡罗样本是从具有这些指定参数值的零假设分布中抽取的。然而,在评估统计量之前,这些参数和零假设分布族的所有其他未知参数始终会拟合到样本,无论是观测 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
是蒙特卡罗样本的数组(具有兼容的形状),axis
是必须计算统计量的data
的轴。- 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 是一个 int,则使用一个新的
RandomState
实例,并使用 random_state 进行种子化。如果 random_state 已经是一个
Generator
或RandomState
实例,则使用该实例。
在 1.15.0 版本中更改:作为从使用
numpy.random.RandomState
到numpy.random.Generator
的 SPEC-007 过渡的一部分,此关键字从 random_state 更改为 rng。在过渡期间,两个关键字将继续有效,尽管一次只能指定一个关键字。在过渡期之后,使用 random_state 关键字的函数调用将发出警告。random_state 和 rng 的行为如上所述,但新代码中应仅使用 rng 关键字。
- dist
- 返回值:
- resGoodnessOfFitResult
具有以下属性的对象。
- fit_result
FitResult
表示提供的 dist 拟合到 data 的对象。此对象包括完全定义零假设分布的分布族参数的值,也就是说,从中抽取蒙特卡罗样本的分布。
- statisticfloat
将提供的 data 与零假设分布进行比较的统计量的值。
- pvaluefloat
零分布中统计量值至少与提供的 data 的统计量值一样极端的元素的比例。
- null_distributionndarray
从零假设分布中抽取的每个蒙特卡罗样本的统计量的值。
- fit_result
注释
这是一个广义的蒙特卡罗拟合优度程序,其中的特殊情况对应于各种 Anderson-Darling 检验、Lilliefors 检验等。该测试在 [2]、[3] 和 [4] 中被描述为参数引导测试。这是一个蒙特卡罗测试,其中指定抽取样本的分布的参数是从数据中估计的。我们在整个过程中使用“蒙特卡罗”而不是“参数引导”来描述该测试,以避免与更熟悉的非参数引导相混淆,并在下面描述如何执行该测试。
传统的拟合优度测试
传统上,对应于一组固定显著性水平的临界值是使用蒙特卡罗方法预先计算的。用户通过仅计算其观测 data 的检验统计量的值,并将此值与列表临界值进行比较来执行该测试。这种做法不是很灵活,因为并非所有分布以及已知和未知参数值的组合都有表格可用。此外,当从有限的列表数据中插值临界值以对应于用户的样本大小和拟合参数值时,结果可能不准确。为了克服这些缺点,此函数允许用户执行适用于其特定数据的蒙特卡罗试验。
算法概述
简而言之,此例程执行以下步骤
将未知参数拟合到给定的 data,从而形成“零假设”分布,并计算此数据和分布对的统计量。
从此零假设分布中抽取随机样本。
将未知参数拟合到每个随机样本。
计算每个样本与其已拟合分布之间的统计量。
将与 data 对应的统计量的值与与 (4) 中随机样本对应的统计量的值进行比较。p 值是统计量值大于或等于观测数据统计量值的样本比例。
更详细地说,步骤如下。
首先,使用最大似然估计将 dist 指定的分布族的任何未知参数拟合到提供的 data。(一个例外是具有未知位置和尺度的正态分布:我们使用偏差校正的标准差
np.std(data, ddof=1)
作为 [1] 中推荐的尺度。)参数的这些值指定了分布族的特定成员,称为“零假设分布”,也就是说,在零假设下抽取数据的分布。计算 data 和零假设分布之间的 statistic,该统计量比较数据与分布。接下来,从零假设分布中抽取许多(具体为 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 中的大多数分布,分布拟合是通过数值优化执行的。
反模式
因此,可能会很想将预先拟合到 data 的分布参数(由用户)视为 known_params,因为指定分布的所有参数排除了将分布拟合到每个蒙特卡罗样本的需要。(这本质上是执行原始 Kilmogorov-Smirnov 检验的方式。)虽然这样的检验可以提供反对零假设的证据,但该检验是保守的,因为小的 p 值往往会(大大)高估犯第一类错误的概率(也就是说,即使零假设为真,也拒绝零假设),并且检验的功效很低(也就是说,即使零假设为假,也不太可能拒绝零假设)。这是因为蒙特卡罗样本不太可能与零假设分布以及 data 一致。这往往会增加零分布中记录的统计量的值,从而使更多的统计量超过 data 的统计量值,从而夸大 p 值。
参考文献
[1] (1,2,3,4,5)M. A. Stephens (1974)。“EDF Statistics for Goodness of Fit and Some Comparisons.” Journal of the American Statistical Association, Vol. 69, pp. 730-737。
[2]W. Stute, W. G. Manteiga, and M. P. Quindimil (1993)。“Bootstrap based goodness-of-fit-tests.” Metrika 40.1: 243-256。
[3]C. Genest, & B Rémillard. (2008)。“Validity of the parametric bootstrap for goodness-of-fit testing in semiparametric models.” Annales de l’IHP Probabilités et statistiques. Vol. 44. No. 6。
[4]I. Kojadinovic and J. Yan (2012)。“Goodness-of-fit testing based on a weighted bootstrap: A fast large-sample alternative to the parametric bootstrap.” Canadian Journal of Statistics 40.3: 480-500。
[5]B. Phipson and G. K. Smyth (2010)。“Permutation P-values Should Never Be Zero: Calculating Exact P-values When Permutations Are Randomly Drawn.” Statistical Applications in Genetics and Molecular Biology 9.1。
[6]H. W. Lilliefors (1967)。“On the Kolmogorov-Smirnov test for normality with mean and variance unknown.” Journal of the American statistical Association 62.318: 399-402。
[7]Filliben, James J. “The probability plot correlation coefficient test for normality.” 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 检验,将观测数据的经验分布函数与正态分布的(理论)累积分布函数进行比较。当然,要做到这一点,必须完全指定零假设下的正态分布。通常的做法是首先将分布的
loc
和scale
参数拟合到观测数据,然后执行检验。>>> 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()
如果分布没有尽可能好地拟合到观测数据,则检验可能无法控制第一类错误率,也就是说,即使零假设为真,也拒绝零假设的可能性。
我们还应该在零分布中寻找可能由不可靠的拟合引起的极端离群值。这些不一定会使结果失效,但它们往往会降低检验的功效。
>>> _, 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()
这个图似乎令人放心。
如果
fit
方法工作可靠,并且如果检验统计量的分布对拟合参数的值不是很敏感,那么预期goodness_of_fit
提供的 p 值是一个很好的近似值。>>> res.statistic, res.pvalue (0.2231991510248692, 0.0525)