非线性回归中的标准误差

Question

我最近在用Python做一些蒙特卡洛物理模拟，但我不知道怎么计算非线性最小二乘拟合中系数的标准误差。

一开始，我以为我的模型是线性的，所以用了SciPy里的scipy.stats.linregress，结果发现其实它是某种幂函数。于是我改用了NumPy的polyfit，设置自由度为2，但我找不到办法来计算系数的标准误差。

我知道gnuplot可以帮我计算误差，但我需要对30多个不同的情况进行拟合。我在想有没有办法让Python从gnuplot读取标准误差，或者有没有其他库可以使用？

Answer 1

看起来gnuplot使用的是一种叫做levenberg-marquardt的方法，而在Python中也有一个实现。你可以通过mpfit.covar这个属性来获取误差估计（顺便提一下，你应该关注这些误差估计“意味着什么”——比如说，其他参数是否可以调整来进行补偿？）

Answer 2

终于找到了这个长期以来被问到的问题的答案！希望这能至少为某些人节省几个小时无望的研究时间。Scipy有一个特别的函数叫做curve_fit，属于它的优化部分。这个函数使用最小二乘法来确定系数，最棒的是，它还会给你一个协方差矩阵。协方差矩阵包含了每个系数的方差。更准确地说，矩阵的对角线就是方差，通过对这些值开平方，就能得到每个系数的标准误差！Scipy对此的文档不多，所以这里有一段示例代码，帮助更好地理解：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plot


def func(x,a,b,c):
    return a*x**2 + b*x + c #Refer [1]

x = np.linspace(0,4,50)
y = func(x,2.6,2,3) + 4*np.random.normal(size=len(x)) #Refer [2] 


coeff, var_matrix = curve_fit(func,x,y)
variance = np.diagonal(var_matrix) #Refer [3]

SE = np.sqrt(variance) #Refer [4]

#======Making a dictionary to print results========
results = {'a':[coeff[0],SE[0]],'b':[coeff[1],SE[1]],'c':[coeff[2],SE[2]]}

print "Coeff\tValue\t\tError"
for v,c in results.iteritems():
    print v,"\t",c[0],"\t",c[1]
#========End Results Printing=================

y2 = func(x,coeff[0],coeff[1],coeff[2]) #Saves the y values for the fitted model

plot.plot(x,y)
plot.plot(x,y2)

plot.show()

1. 这个函数返回的内容非常重要，因为它定义了将用于拟合模型的内容
2. 使用这个函数创建一些任意的数据加上一些噪声
3. 将协方差矩阵的对角线保存到一个一维矩阵中，这其实就是一个普通的数组
4. 对方差开平方以获得标准误差（SE）