如何使用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检验 是非常优雅的。

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

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

在statsmodel中有一个 卡方拟合优度函数，但它假设数据是离散分布（而指数分布是连续的）。

官方的 scipy.stats教程 只涵盖了自定义分布的情况，概率是通过调整很多表达式（npoints, npointsh, nbound, normbound）来构建的，所以我不太清楚如何处理其他分布。卡方示例假设期望值和自由度已经得到了。

另外，我并不是想要“手动”执行检验，正如 这里讨论过的那样，我想知道如何使用现有的库函数。

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 格式化，并调整成1维的10,000个点。

现在我们来用卡方检验（ChiSquared test）来评估这个拟合效果：

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

你为什么需要“验证”它是否是指数型的呢？你确定需要做统计测试吗？我几乎可以保证它最终不是指数型的，如果你有足够的数据，测试结果会很显著，这样一来，使用这个测试的逻辑就显得有些勉强了。你可以看看这个讨论：正态性测试“基本上没用”吗？，或者我在这里的回答：用大量观察数据测试异方差性。

通常来说，使用qq图和/或pp图会更好（具体取决于你是关心分布的尾部还是中间部分的拟合情况，看看我在这里的回答：PP图与QQ图的比较）。关于如何在Python的SciPy中制作qq图的信息，可以在这个讨论中找到：使用SciPy绘制分位数-分位数图

Answer 3

这里有一个关于如何让每个区间的概率差不多的简单方法：

• 首先，估算一下数据的分布参数。
• 接着，如果你用的是scipy.stats这个库，可以用反累积分布函数（inverse cdf），也叫ppf，来得到一个均匀概率网格的区间边界。比如，你可以这样写：distribution.ppf(np.linspace(0, 1, n_bins + 1), *args)。
• 然后，使用np.histogram来统计每个区间里有多少个数据。

接下来，可以对这些频率使用卡方检验。

另外一种方法是从排序后的数据中找到百分位数，然后用累积分布函数（cdf）来找出实际的概率。

需要注意的是，这种方法只是个近似，因为卡方检验的理论是基于在分组数据上用最大似然法估算参数的。我不太确定根据数据选择区间边界是否会影响渐近分布。

我很久没研究这个了。如果这个近似方法不够好，建议你去stats.stackexchange上提问。