我有一组实验值和一个概率密度函数,它可以描述它们的分布:

def bekkers(x, a, m, d):
    p = a*np.exp((-1*(x**(1/3) - m)**2)/(2*d**2))*x**(-2/3)
    return(p)

我用scipy.optimize.curve_fit估计了函数的参数,现在我需要测试拟合优度。我发现了一个scipy.stats.kstest函数,它可以很好地满足我的需求,但它需要一个连续分布函数。我如何完成我的任务?

最佳答案

注:我不确定你可能的x值的范围是什么,以及你对a,m,d的估计是什么,所以我尽量让这些尽可能开放。
KS测试的CDF指的是一个累积分布,而不是一个连续分布函数(您已经得到了)。我们将为此建立一个函数,因为我不确定你所提供的方程的积分是否有一个封闭形式,所以我们将用scipy.integrate来实现。
为了与其他numpy/scipy工具一起使用,我们希望它接收并返回一个数组(可能有一种更漂亮的方法可以这样做,但下面的方法仍然有效)。还要注意,必须规范化cdf,因为至少对于我选择的值和范围,可能值的整个范围内的积分不等于1。下面是它的外观:

def bekkers_cdf(x,a,m,d,range_start,range_end):
    values = []
    for value in x:
        integral = integrate.quad(lambda k: bekkers(k,a,m,d),range_start,value)[0]
        normalized = integral/integrate.quad(lambda k: bekkers(k,a,m,d),range_start,range_end)[0]
        values.append(normalized)
    return np.array(values)

一旦我们有了这个,我们现在可以评估我们的ks.测试(使用我为范围和a,m,d所做的一些值):
my_start,my_end = 1,10
my_a,my_m,my_d = 1,1,1
my_data = [1.5,1.6,1.8,2.1,2.2,3.3,4,6,8,9]
stats.kstest(my_data,lambda x: bekkers_cdf(x,my_a,my_m,my_d,my_start,my_end))

这将返回:
(0.17609125905568074, 0.9157727421346824)

第一个值是统计值,第二个是p值。有这么高的p值,我们绝对不能拒绝这个数据来自这个分布。
代码摘要:
import numpy as np
import scipy as sp
from scipy import integrate,stats

def bekkers(x, a, m, d):
    p = a*np.exp((-1*(x**(1/3) - m)**2)/(2*d**2))*x**(-2/3)
    return(p)

def bekkers_cdf(x,a,m,d,range_start,range_end):
    values = []
    for value in x:
        integral = integrate.quad(lambda k: bekkers(k,a,m,d),range_start,value)[0]
        normalized = integral/integrate.quad(lambda k: bekkers(k,a,m,d),range_start,range_end)[0]
        values.append(normalized)
    return np.array(values)

my_start = 1
my_end = 10
my_a,my_m,my_d = 1,1,1
my_data = [1.5,1.6,1.8,2.1,2.2,3.3,4,6,8,9]
stats.kstest(my_data,lambda x: bekkers_cdf(x,my_a,my_m,my_d,my_start,my_end))

为了找点乐子,我们可以看一下ks测试在看什么。为了做到这一点,我们绘制了理论cdf从我们的数据相比,提出的函数。(注意,下面我在数据的cdf点中硬编码,但这很容易编程)使用matplotlib这是:
import matplotlib.pyplot as plt
xs = np.linspace(1, 10)
ys = bekkers_cdf(xs,my_a,my_m,my_d,my_start,my_end)
theoretical, =plt.plot(xs,ys,linewidth=2)
x2s = [1,1.5,1.6,1.8,2.1,2.2,3.3,4,6,8,9,10]
y2s = [0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1,1]
data, =plt.plot(x2s,y2s,linewidth=2)
plt.legend([theoretical,data],['theoretical','data'])

结果是:
我们看到数据的cdf与建议的分布相似,因此我们的测试没有拒绝来自此分布的样本数据的空值是有意义的。

关于python - 如何在python中的自定义概率密度函数上执行Kolmogorov-Smirnov拟合优度检验?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/26143499/

10-09 08:32
查看更多