如何在Python中计算最小二乘拟合(scipy.optimize.leastsq)的置信区间?
如何在Python中计算最小二乘拟合(scipy.optimize.leastsq)的置信区间?
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
我不确定您所说的置信区间是什么意思。
一般来说,leastsq
并不知道您试图最小化的函数很多细节,因此它无法真正给出置信区间。然而,它确实返回了Hessian的估计值,也就是将二阶导数推广到多维问题的方法。
如函数的docstring中所示,您可以使用该信息以及残差(拟合解与实际数据之间的差异)来计算参数估计的协方差,这是置信区间的局部猜测。
请注意,这只是一个局部信息,我怀疑您只有在目标函数严格凸时才能做出严格的结论。我没有任何证明或参考资料支持这个说法 :)