python - 如何在 python 中对自定义概率密度函数执行 Kolmogorov-Smirnov 拟合优度检验？

Question

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

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 函数，它可以完全满足我的需要，但它需要一个连续的分布函数。我如何处理我的任务？

Answer 1

注意:我不确定你可能的 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.test(使用我为范围、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-test 在看什么。为此，我们将我们数据的理论 cdf 与建议函数的 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 与建议分布所建议的 cdf 相似，因此我们的测试没有拒绝 null 且样本数据来自该分布是有道理的。