18

假设我有一些经验数据：

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 中有一个卡方 GoF 函数，但它假设离散分布（并且指数分布是连续的）。

官方 scipy.stats 教程仅涵盖自定义分布的情况，并且概率是通过摆弄许多表达式（npoints、npointsh、nbound、normbound）来构建的，所以我不太清楚如何为其他分布做这件事。卡方示例假设已经获得了预期值和自由度。

另外，我不是在寻找一种“手动”执行测试的方法，正如这里已经讨论过的那样，而是想知道如何应用一个可用的库函数。

4

3 回答 3

4

等概率箱的近似解：

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

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

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

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

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

于 2014-06-24T14:05:28.820 回答
2

为什么你需要“验证”它是指数级的？你确定你需要统计测试吗？我几乎可以保证这最终不是指数级的，如果你有足够的数据，测试会很重要，这使得使用测试的逻辑相当强制。阅读此 CV 主题可能会对您有所帮助：正态性测试“基本上没用”吗？，或者我在这里的回答：用许多观察结果检验异方差性。

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

于 2014-06-24T14:25:48.993 回答
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)

在此处输入图像描述

于 2020-11-05T14:54:34.490 回答