在scipy python中用leastsq拟合的置信区间

Question

如何在Python中计算最小二乘拟合的置信区间（使用scipy.optimize.leastsq）？

Answer 1

估算置信区间（CI）最简单的方法就是把标准误差（标准差）乘以一个常数。要计算这个常数，你需要知道自由度（DOF）和你想要计算的置信水平。通过这种方法估算的置信区间，有时被称为渐近置信区间。

你可以在Motulsky和Christopoulos的书《使用线性和非线性回归拟合生物数据》中了解更多内容，书的链接在这里：google books。同样的书（或者非常相似的书）可以在这里免费获取：作为作者软件的手册。

你还可以阅读一下如何使用C++ Boost.Math库计算置信区间。在这个例子中，置信区间是针对一个变量的分布计算的。在最小二乘拟合的情况下，自由度不是-1，而是，其中是参数的数量。在Python中做同样的事情应该也很简单。

这就是最简单的估算方法。我不太了解zephyr提出的自助法，但它可能比我写的这种方法更可靠。

Answer 2

我不太明白你说的置信区间是什么意思。

一般来说，leastsq这个函数对你想要最小化的那个函数了解得不多，所以它其实不能给出置信区间。不过，它会返回一个关于Hessian的估计，简单来说就是把二阶导数的概念扩展到多维问题上。

正如函数的文档中提到的，你可以利用这些信息和残差（也就是你拟合的结果和实际数据之间的差异）来计算参数估计的协方差，这可以算是对置信区间的一个局部估计。

需要注意的是，这只是局部的信息，我怀疑只有在你的目标函数是严格凸的情况下，才能严格得出结论。我对此没有任何证明或参考资料 :）。

Answer 3

我会使用自助法（bootstrapping）。
可以看看这里: 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

在scipy python中用leastsq拟合的置信区间

3 个回答

撰写回答