我有一个实验数据数据集(y_observed,t,其中 y_observed 表示测量值,t 表示自测量开始以来的时间(以秒为单位))。我对这些数据进行高斯拟合,但想要评估拟合高斯是否有意义或者观察到的实验数据是否是非高斯的。
我尝试在 Python 中实现 2 个样本 Kolmogorov-Smirnov 拟合优度检验。然而,我想知道我是否正确地实施和解释了这个测试。下面是一个小例子,其中 y_observed 看起来不像高斯分布(见附图),但根据 2 个样本 K-S 检验 y_observed 确实遵循高斯分布。有人可以向我解释我做错了什么吗?预先感谢!
from scipy.stats import ks_2samp
import numpy as np
y_observed = [
1.93291, 9.84869, 29.5713, 20.6346, -1.51537,
0.396292, 50.8895, 68.855, 63.9291, 55.6863,
37.6503, 46.87, 33.6637, 25.0395,18.3027,
38.0947, 27.9305, 10.2012, -0.704519, 18.6656,
20.2873, 8.78955, 4.39672
]
t = [
519, 522, 525, 528, 531,
534, 537, 540, 543, 546,
549, 552, 555, 558 ,561,
564, 567, 570, 573, 576,
579,582,585
]
y_fit = [
5.60417, 8.74951, 12.9907, 18.3426, 24.63,
31.4518, 38.1947, 44.1102, 48.4454, 50.5991,
50.2586, 47.4739, 42.6458, 36.4314, 29.5973,
22.8668, 16.8011, 11.7394, 7.80065, 4.92939,
2.96233, 1.69298, 0.920122
]
count, bins_count = np.histogram(y_observed, bins=5)
pdf = count / sum(count)
cdf = np.cumsum(pdf)
count_2, bins_count_2 = np.histogram(y_fit, bins=5)
pdf_2 = count_2 / sum(count_2)
cdf_2 = np.cumsum(pdf_2)
KS_stat, P_val = ks_2samp(cdf, cdf_2)
n_1, n_2 = len(y_observed), len(y_fit)
critical_val = np.sqrt((n_1+n_2)/(n_1*n_2))*1.36
if P_val<0.05 and KS_stat>critical_val:
print("Observed data does not follow a Gaussian distribution")
else:
print("Observed data follows a Gaussian distribution")
有人可以向我解释一下我做错了什么吗?
stats.ks_2samp
接受两个样本,并对这两个样本来自相同(未知)连续分布的原假设进行检验。因此,第一个问题是代码没有传递两个样本 - 它似乎传递了两个数组的粗粒度经验 CDF,其中一个是您的样本。 (我不知道 y_fit
是。)如果你想评估 y_observed
和 y_fit
是否来自同一分布(如果你想知道 y_observed
是否被绘制,这看起来没有意义)从正态分布),你只需调用
ks_2samp(y_observed, y_fit)
第二个问题是,这不是评估数据是否非高斯的适当测试。
如果必须进行 KS 测试,您需要单样本 KS 测试。
from scipy import stats
stats.ks_1samp(y_observed, stats.norm().cdf)
问题是,正如所写的,这将测试数据是从“标准”正态分布中提取的零假设。如果你知道假设正态分布的参数loc
和
scale
,你可以将它们传入:from scipy import stats
stats.ks_1samp(y_observed, stats.norm(loc=loc, scale=scale).cdf)
如果您不知道参数(即使您可以将它们拟合到数据),那么此测试是不合适的。
SciPy 中还有其他几个您可以使用的正态性测试。您只需提供要测试的样本即可通过它们。
stats.shapiro(y_observed)
stats.jarque_bera(y_observed)
stats.normaltest(y_observed)
stats.kurtosistest(y_observed)
stats.skewtest(y_observed)
stats.anderson(y_observed)
通常是首选。
anderson
将是另一个不错的选择,但是 scipy.stats.anderson
不会产生 p 值(它只返回临界值表),因此它不太容易使用。 scipy.stats.goodness_of_fit
是最灵活的(并且可以为 Anderson-Darling 检验生成 p 值),但您不需要它来进行正态性检验。请注意,无论您选择什么,频率主义零假设检验都不可能得出数据是从高斯分布中提取的结论。它只能提供证据表明数据不是从高斯分布中得出的;否则是不确定的。