从scipy.optimize.curve_fit获取参数估计的标准误差

Question

我正在使用 scipy.optimize.curve_fit 来给一些数据拟合曲线。大部分情况下，拟合效果看起来不错。不过，当我打印出结果时，发现 pcov 的值是无穷大（inf）。

我其实需要计算与我拟合的参数相关的误差，但即使它给了我协方差矩阵，我也不太确定该怎么做。

我拟合的模型是：

def intensity(x,R_out,R_in,K_in,K_out,a,b,c):
    K_in,K_out = abs(0.0),abs(K_out)
    if x<=R_in:
        return 2*R_out*(K_out*np.sqrt(1-x**2/R_out**2)-
                (K_out-0.0)*np.sqrt(R_in**2/R_out**2-x**2/R_out**2)) + c
    elif x>=R_in and x<=R_out:
        return K_out*2*R_out*np.sqrt(1-x**2/R_out**2) + c
    elif x>R_out:
        return c

intensity_vec = np.vectorize(intensity)



def intensity_vec_self(x,R_out,R_in,K_in,K_out,a,b,c):
    y = np.zeros(x.shape)
    for i in range(len(y)):
        y[i]=intensity_vec(x[i],R_out,R_in,K_in,K_out,a,b,c)
    return y

我有400个数据点，如果你觉得有帮助，我可以把这些数据放上来。

总之，我无法让 curve_fit 打印出我的 pcov，需要帮助来弄清楚为什么会这样，以及我是否能让它正常工作。

另外，如果能简单解释一下，我也想知道如何使用 pcov 数组来获得与我的拟合相关的误差。

谢谢！

Answer 1

参数的方差就是方差-协方差矩阵中的对角线元素，而标准误差则是这些方差的平方根。你可以用这个代码来计算：np.sqrt(np.diag(pcov))

关于出现inf的情况，可以看看这两个例子进行比较：

In [129]:
import numpy as np
def func(x, a, b, c, d):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5, 1)
ydata = y + 0.2 * np.random.normal(size=len(xdata))
popt, pcov = so.curve_fit(func, xdata, ydata)
print np.sqrt(np.diag(pcov))
[ inf  inf  inf  inf]

还有：

In [130]:

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
ydata = y + 0.2 * np.random.normal(size=len(xdata))
popt, pcov = so.curve_fit(func, xdata, ydata)
print np.sqrt(np.diag(pcov))
[ 0.11097646  0.11849107  0.05230711]

在这个极端的例子中，d对函数func没有任何影响，因此它的方差会是+inf，换句话说，它可以是任何值。把d从func中去掉，结果就会变得合理。

实际上，如果参数的量级差别很大，比如：

def func(x, a, b, c, d):
    #return a * np.exp(-b * x) + c
    return a * np.exp(-b * x) + c + d*1e-10

你也会因为浮点数溢出或下溢而得到inf。

在你的情况中，我觉得你从来没有使用过a和b。这就像第一个例子一样。