如何使用 Python 中的科学库执行卡方拟合优度测试？

Question

假设我有一些经验数据：

from scipy import stats
size = 10000
x = 10 * stats.expon.rvs(size=size) + 0.2 * np.random.uniform(size=size)

它呈指数分布（带有一些噪声），我想使用卡方拟合优度 (GoF) 测试来验证这一点。使用 Python 中的标准科学库（例如 scipy 或 statsmodels）以最少的手动步骤和假设来执行此操作的最简单方法是什么？

我可以为模型拟合：

param = stats.expon.fit(x)
plt.hist(x, normed=True, color='white', hatch='/')
plt.plot(grid, distr.pdf(np.linspace(0, 100, 10000), *param))

计算Kolmogorov-Smirnov test非常优雅。

>>> stats.kstest(x, lambda x : stats.expon.cdf(x, *param))
(0.0061000000000000004, 0.85077099515985011)

但是，我找不到计算卡方检验的好方法。

有一个chi-squared GoF function in statsmodel，但它假设一个离散分布（并且指数分布是连续的）。

official scipy.stats tutorial 仅涵盖自定义分布的情况，概率是通过摆弄许多表达式（npoints、npointsh、nbound、normbound）来构建的，所以我不太清楚如何为其他分布做这件事。 chisquare examples 假设已获得预期值和景深。

另外，我不是像already discussed here 那样寻找“手动”执行测试的方法，而是想知道如何应用其中一个可用的库函数。

Answer 1

我用 OpenTURNS 试过你的问题。 开头是一样的：

import numpy as np
from scipy import stats
size = 10000
x = 10 * stats.expon.rvs(size=size) + 0.2 * np.random.uniform(size=size)

如果您怀疑您的样本 x 来自指数分布，您可以使用 ot.ExponentialFactory() 来拟合参数：

import openturns as ot
sample = ot.Sample([[p] for p in x])
distribution = ot.ExponentialFactory().build(sample)

由于 Factory 需要一个 ot.Sample() 作为输入，我需要格式化 x 并将其重塑为 10.000 个维度为 1 的点。

现在让我们使用卡方检验评估此拟合：

result = ot.FittingTest.ChiSquared(sample, distribution, 0.01)
print('Exponential?', result.getBinaryQualityMeasure(), ', P-value=', result.getPValue())
>>> Exponential? True , P-value= 0.9275212544642293

很好！

当然，print(distribution) 会为您提供拟合参数：

>>> Exponential(lambda = 0.0982391, gamma = 0.0274607)

Answer 2

为什么需要“验证”它是指数级的？你确定你需要统计测试吗？我几乎可以保证这最终不是指数级的，如果你有足够的数据，测试会很重要，这使得使用测试的逻辑相当强制。阅读此 CV 线程：Is normality testing 'essentially useless'?，或我的答案：Testing for heteroscedasticity with many observations，可能会对您有所帮助。

通常最好使用 qq-plot 和/或 pp-plot（取决于您是否关心分布的尾部或中间的拟合，请参阅我的答案：PP-plots vs. QQ-plots）。关于如何在 Python SciPy 中制作 qq-plots 的信息可以在这个 SO 线程中找到：Quantile-Quantile plot using SciPy

Answer 3

等概率箱的近似解：

• 估计分布参数
• 如果是 scipy.stats.distribution，则使用逆 cdf ppf 来获取规则概率网格的边，例如distribution.ppf(np.linspace(0, 1, n_bins + 1), *args)
• 然后，使用 np.histogram 计算每个 bin 中的观察次数

然后对频率使用卡方检验。

另一种方法是从已排序数据的百分位数中查找 bin 边缘，并使用 cdf 查找实际概率。

这只是近似值，因为卡方检验的理论假设参数是通过分箱数据的最大似然估计的。而且我不确定基于数据选择binedges是否会影响渐近分布。

我很久没有研究这个问题了。 如果一个近似的解决方案不够好，那么我建议你在 stats.stackexchange 上提问。