在Python的Scipy库中,使用leastsq拟合方法计算置信区间。

7

如何在Python中计算最小二乘拟合(scipy.optimize.leastsq)的置信区间?


你可以使用Bootstrap:https://dev59.com/questions/_6Lia4cB1Zd3GeqPh1to#66008548 - Marco Cerliani
3个回答

9
我会使用引导法。
请参考这里:http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html 噪声高斯的简单示例:
x = arange(-10, 10, 0.01)

# model function
def f(p):
    mu, s = p
    return exp(-(x-mu)**2/(2*s**2))

# create error function for dataset    
def fff(d):
    def ff(p):
        return d-f(p)
    return ff

# create noisy dataset from model
def noisy_data(p):
    return f(p)+normal(0,0.1,len(x))

# fit dataset to model with least squares    
def fit(d):
    ff = fff(d)
    p = leastsq(ff,[0,1])[0]
    return p

# bootstrap estimation        
def bootstrap(d):
    p0 = fit(d)
    residuals = f(p0)-d
    s_residuals = std(residuals)

    ps = []
    for i in range(1000):
        new_d = d+normal(0,s_residuals,len(d))
        ps.append(fit(new_d))

    ps = array(ps)
    mean_params = mean(ps,0)
    std_params = std(ps,0)

    return mean_params, std_params

data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)

print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996

2
这是一个非常好的答案,如果您能够包括几个简单的英语句子来解释什么是引导程序以及它如何工作,那么它将受益匪浅。 - jb.

4

我不确定您所说的置信区间是什么意思。

一般来说,leastsq并不知道您试图最小化的函数很多细节,因此它无法真正给出置信区间。然而,它确实返回了Hessian的估计值,也就是将二阶导数推广到多维问题的方法。

如函数的docstring中所示,您可以使用该信息以及残差(拟合解与实际数据之间的差异)来计算参数估计的协方差,这是置信区间的局部猜测。

请注意,这只是一个局部信息,我怀疑您只有在目标函数严格凸时才能做出严格的结论。我没有任何证明或参考资料支持这个说法 :)


2
估计置信区间(CI)最简单的方法是将标准误差(标准偏差)乘以一个常数。要计算常数,您需要知道自由度(DOF)的数量和要计算CI的置信水平。 用这种方式估计的CI有时被称为渐近CI。 您可以在Motulsky&Christopoulos的“使用线性和非线性回归将模型拟合到生物数据”中阅读更多内容(谷歌图书)。该书(或非常相似的书籍)可免费作为作者软件的手册获得。
您还可以阅读如何使用C ++ Boost.Math库计算CI。在此示例中,CI是针对一个变量的分布计算的。在最小二乘拟合的情况下,DOF不是N-1,而是N-M,其中M是参数的数量。在Python中完成相同的操作应该很容易。
这是最简单的估计方法。我不知道zephyr提出的自举方法,但它可能比我所写的方法更可靠。

问题是如何计算它们,以便使用Python函数访问它们? - rhody

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接