Python单面KS-Test

Question

我已经针对大量理论概率分布对我的分布进行了单边 KS 测试(观察公共(public)交通公共(public)交通网络的占用情况，值从 0 到 100):

cdfs = [
    "norm",            #Normal (Gaussian)
    "alpha",           #Alpha
    "anglit",          #Anglit
    "arcsine",         #Arcsine
    "beta",            #Beta
    "betaprime",       #Beta Prime
    "bradford",        #Bradford
    "burr",            #Burr
    "cauchy",          #Cauchy
     ....
]

for cdf in cdfs:
    #fit our data set against every probability distribution
    parameters = eval("scipy.stats."+cdf+".fit(data_sample)");

    #Applying the Kolmogorov-Smirnof one sided test
    D, p = scipy.stats.kstest(data_sample, cdf, args=parameters);

    #pretty-print the results
    print (cdf.ljust(16) + ("p: "+str('{0:.10f}'.format(p)).ljust(40)+"D: "+str('{0:.10f}'.format(D))));

根据我对单面 KS-Test 的理解，最适合我的数据的理论分布是单面 KS-Test 返回大 p 值和低 D-KSstatistic 值的分布。

据此，最佳拟合是:

cdf: invweibull p:0.1624542096  D:0.0352622822
cdf: genextreme p:0.1624292228  D:0.0352633673
cdf: nct        p:0.1280588168  D:0.0369024688
cdf: invgamma   p:0.1273446642  D:0.0369401507
cdf: johnsonsu  p:0.0449026953  D:0.0433976894
cdf: invgauss   p:0.0336248605  D:0.0450259762
(...)
cdf: frechet_l  p:0.0000000000  D:0.8405035144
cdf: reciprocal p:0.0000000000  D:0.9380000000
cdf: truncnorm  p:0.0000000000  D:0.9380000000
cdf: powernorm  p:0.0000000000  D:1.0000000000

此外，当我尝试将这些所谓的最佳拟合分布直观地拟合到我的数据时，有些东西并没有加起来:

from scipy.stats import invgauss, invweibull, genextreme

fig, ax = plt.subplots(1, 1)
mu = 10.145462645553
x = np.linspace(invgauss.ppf(0.75, mu), invgauss.ppf(0.975, mu), 100)
ax.plot(x, invgauss.pdf(x, mu), 'r-', color='green', lw=1, alpha=0.6, label='invgauss pdf')

c = 0.8
y = np.linspace(invweibull.ppf(0.75, c), invweibull.ppf(0.975, c), 100)
ax.plot(y, invweibull.pdf(y, c), 'r-', color='red', lw=1, alpha=0.6, label='invweibull pdf')

c = -1.5
z = np.linspace(genextreme.ppf(0.75, c), genextreme.ppf(0.96, c), 100)
ax.plot(z, genextreme.pdf(z, c), 'r-', lw=1, color='yellow', alpha=0.6, label='genextreme pdf')

ax.hist(data_sample, normed=True, histtype='stepfilled', bins=20, alpha=0.2, label='my distribution')
ax.legend(loc='best', frameon=False)
plt.show()

结果似乎与我的数据不匹配 invgauss、invweibull 或 genextreme 概率分布。



我做错了什么或假设 KS 测试结果有问题吗？

我的分布中的数据样本:

array([ 29.75,   0.8 ,   9.  ,   4.77,  28.75,  31.1 ,  52.12,   5.  ,
        10.55,  17.26,  19.28,  25.77,  53.13,  28.  ,   4.1 ,   2.92,
        40.4 ,  15.33,  10.62,  20.6 ,  26.11,  15.  ,   5.3 ,  38.87,
         1.28,   1.5 ,  20.88,  16.  ,  10.33,   6.5 ,   6.  ,  22.5 ,
         7.88,   2.72,  60.33,  26.14,  18.  ,  18.58,  25.  ,  69.62,
         0.5 ,   0.  ,  26.87,  11.85,  13.16,  39.45,  17.6 ,  14.66,
        84.52,   3.62,  30.33,   4.25,  25.  ,  35.  ,  28.85,  48.37,
        12.55,  50.  ,  22.94,   7.42,   2.37,  49.66,  22.94,   7.57,
       101.12,   4.42,  43.88,   7.  ,  13.  ,  31.12,  20.71,   0.  ,
        22.  ,  21.34,  23.61,   0.5 ,  16.23,  27.11,   2.22,  59.  ,
        24.41,  41.69,   2.68,  49.  ,  51.6 ,  95.8 ,   0.  ,  26.8 ,
        66.  ,  43.02,  13.85,  46.91,  38.77,   6.5 ,  24.  ,  54.14,
        50.81,  21.55,  19.22,  12.83])

解决方案

请参阅接受的答案以获取更多详细信息。仅供引用，在估计正确的参数并将其传递给单边 KS 测试认为与我自己的分布相似的最相似的理论分布后，我能够直观地确认分布相似性。

Answer 1

简答

你说清楚了，只留下一件事:不同的分布有不同的参数 .我们应该将估计的参数传递给分布，然后执行 KS-test 和您的最终密度图。

scipy.stats.invgamma.fit(data_sample),\
scipy.stats.norm.fit(data_sample)

((4.399779777260058, -15.382411650381744, 137.60256212682822),
 (24.501099999999997, 21.016423572768037))

换句话说，如果你想用各种分布测试你的数据，你应该仔细为每个分布设置参数。

首先，您将数据与分布拟合并获得每个分布的估计参数。

接下来，您对估计的分布执行 KS 检验(在第一步中使用拟合参数)。

最后，你应该绘制估计的分布(应该将参数传递给每个分布)和你的原始数据，看看 KS-test 的结果是否可信。

修订代码

from scipy.stats import bradford,invgauss, invweibull, genextreme

fig, ax = plt.subplots(1, 1)

# set colors for different distributions
colors = ['purple','green','red','yellow'] 
# the distributions you want to compare, add a braford to see if there is any difference
dists = [bradford,invgauss, invweibull, genextreme] 
for cdf in dists:
    # get the paras. Note that the order of parameters is the same with .ppf(parameters), due to the consitiancy of scipy.stats
    parameters =  eval("scipy.stats."+cdf.name+".fit(data_sample)") 
    x = np.linspace(cdf.ppf(0.3, *parameters), cdf.ppf(0.99, *parameters), 100)
    ax.plot(x, cdf.pdf(x, *parameters), 'r-', color=colors[dists.index(cdf)], lw=1, alpha=0.6, label= cdf.name + ' pdf')

ax.hist(data_sample, density=True, histtype='stepfilled', bins=20, alpha=0.2, label='my distribution')
ax.legend(loc='best', frameon=False)
plt.show()

笔记

由于我们使用 scipy.stats用于参数的拟合和 ppf 生成。保持每个分布中参数的顺序。如果你使用不同的包或软件做同样的事情，参数的顺序不保证是一样的!

例如，在 R 中， Gamma 分布只有两个参数:位置和形状，这与我的统计教科书相同。然而，在 python 中， Gamma 分布有 3 个参数。过去，我必须在 R 中编写自己的 gamma 函数来生成具有 3 个参数的 gamma 分布。

代码:

scipy.stats.gamma.fit(data_sample)

输出:

(0.9205943551269966, -6.9329968733926e-26, 27.77113522169767)

正如卢卡斯在评论中所说， scipy.stats只提供 .fit连续分布的方法。如果要拟合离散分布，请参阅 statsmodels.discrete.discrete_model .职位分配可以使用 MLE 或 矩估计器用于拟合 lambda，例如 lambda = 1/sample_mean。估算方法供您选择。您可以在特定情况下编写自己的方法