Python 单侧 KS 检验
Python one-sided KS-Test
我已经 运行 对我的分布进行了单边 KS 检验（观察质量 t运行sit public t运行sportation grid值 运行 从 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 检验的理解，最适合我的数据的理论分布是单向 KS 检验 returns 大 p 值和低 D 的分布-KS统计值。
据此，最合适的是：
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 测试认为与我自己的分布相似的最相似的理论分布后，我能够直观地确认分布相似性。
简答
你说清楚了，只剩下一件事：不同的分布有不同的参数。我们应该将估计的参数传递到分布中，然后执行 KS 检验和最终的密度图。
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分布只有2个参数：location和shape，这和我的统计学教科书一样。但是，在 python 中，伽玛分布有 3 个参数。过去，我必须在 R 中编写自己的 gamma 函数来生成具有 3 个参数的 gamma 分布。
代码：
scipy.stats.gamma.fit(data_sample)
输出：
(0.9205943551269966, -6.9329968733926e-26, 27.77113522169767)
正如Lucas在评论中所说，scipy.stats只提供了.fit连续分布的方法。如果您想要拟合离散分布，请参阅 statsmodels.discrete.discrete_model。至于possion分布，你可以使用MLE或moment estimator来拟合lambda，例如lambda = 1/sample_mean。估算方法供您选择。您可以在特定情况下编写自己的方法
我已经 运行 对我的分布进行了单边 KS 检验（观察质量 t运行sit public t运行sportation grid值 运行 从 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 检验的理解，最适合我的数据的理论分布是单向 KS 检验 returns 大 p 值和低 D 的分布-KS统计值。
据此，最合适的是：
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 测试认为与我自己的分布相似的最相似的理论分布后，我能够直观地确认分布相似性。
简答
你说清楚了，只剩下一件事：不同的分布有不同的参数。我们应该将估计的参数传递到分布中，然后执行 KS 检验和最终的密度图。
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分布只有2个参数：location和shape，这和我的统计学教科书一样。但是，在 python 中，伽玛分布有 3 个参数。过去，我必须在 R 中编写自己的 gamma 函数来生成具有 3 个参数的 gamma 分布。
代码：
scipy.stats.gamma.fit(data_sample)
输出：
(0.9205943551269966, -6.9329968733926e-26, 27.77113522169767)
正如Lucas在评论中所说，scipy.stats只提供了.fit连续分布的方法。如果您想要拟合离散分布，请参阅 statsmodels.discrete.discrete_model。至于possion分布，你可以使用MLE或moment estimator来拟合lambda，例如lambda = 1/sample_mean。估算方法供您选择。您可以在特定情况下编写自己的方法