在Python中使用optimize.leastsq获取拟合参数的标准误差

Question

我有一组数据（位移和时间的关系），我用 optimize.leastsq 方法对这些数据进行了拟合。现在我想要得到拟合参数的误差值。在查看文档时，我发现输出的矩阵是雅可比矩阵，我需要把这个矩阵和残差矩阵相乘才能得到我的值。不过，我对统计学不是很了解，所以在这些术语中有点迷失。

根据我的理解，我只需要与我的拟合参数相关的协方差矩阵，这样我就可以对对角线上的元素开平方，得到拟合参数的标准误差。我模糊记得协方差矩阵是 optimize.leastsq 方法输出的内容之一。这是对的吗？如果不是，我该如何获取残差矩阵，以便将输出的雅可比矩阵相乘，从而得到我的协方差矩阵呢？

任何帮助都将不胜感激。我对 Python 非常陌生，因此如果这个问题很基础，我在此先道歉。

拟合的代码如下：

fitfunc = lambda p, t: p[0]+p[1]*np.log(t-p[2])+ p[3]*t # Target function'

errfunc = lambda p, t, y: (fitfunc(p, t) - y)# Distance to the target function

p0 = [ 1,1,1,1] # Initial guess for the parameters


  out = optimize.leastsq(errfunc, p0[:], args=(t, disp,), full_output=1)

这里的 t 和 disp 是时间和位移值的数组（基本上就是两列数据）。我在代码的开头导入了所有需要的内容。拟合的值和输出提供的矩阵如下：

[  7.53847074e-07   1.84931494e-08   3.25102795e+01  -3.28882437e-11]

[[  3.29326356e-01  -7.43957919e-02   8.02246944e+07   2.64522183e-04]
 [ -7.43957919e-02   1.70872763e-02  -1.76477289e+07  -6.35825520e-05]
 [  8.02246944e+07  -1.76477289e+07   2.51023348e+16   5.87705672e+04]
 [  2.64522183e-04  -6.35825520e-05   5.87705672e+04   2.70249488e-07]]

我怀疑目前的拟合结果可能有点问题。等我能得到误差值时就能确认这一点。

Answer 1

在进行线性回归时，我们可以准确计算出误差。实际上，leastsq这个函数会给出不同的值：

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt

A = 1.353
B = 2.145

yerr = 0.25

xs = np.linspace( 2, 8, 1448 )
ys = A * xs + B + np.random.normal( 0, yerr, len( xs ) )

def linearRegression( _xs, _ys ):
    if _xs.shape[0] != _ys.shape[0]:
        raise ValueError( 'xs and ys must be of the same length' )
    xSum = ySum = xxSum = yySum = xySum = 0.0
    numPts = 0

    for i in range( len( _xs ) ):
        xSum += _xs[i]
        ySum += _ys[i]
        xxSum += _xs[i] * _xs[i]
        yySum += _ys[i] * _ys[i]
        xySum += _xs[i] * _ys[i]
        numPts += 1

    k = ( xySum - xSum * ySum / numPts ) / ( xxSum - xSum * xSum / numPts )
    b = ( ySum - k * xSum ) / numPts
    sk = np.sqrt( ( ( yySum - ySum * ySum / numPts ) / ( xxSum - xSum * xSum / numPts ) - k**2 ) / numPts )
    sb = np.sqrt( ( yySum - ySum * ySum / numPts ) - k**2 * ( xxSum - xSum * xSum / numPts ) ) / numPts

    return [ k, b, sk, sb ]

def linearRegressionSP( _xs, _ys ):
    defPars = [ 0, 0 ]
    pars, pcov, infodict, errmsg, success = \
    leastsq( lambda _pars, x, y: y - ( _pars[0] * x + _pars[1] ), defPars, args = ( _xs, _ys ), full_output=1 )

    errs = []
    if pcov is not None:
        if( len( _xs ) > len(defPars) ) and pcov is not None:
            s_sq = ( ( ys - ( pars[0] * _xs + pars[1] ) )**2 ).sum() / ( len( _xs ) - len( defPars ) )
            pcov *= s_sq

        for j in range( len( pcov ) ):
            errs.append( pcov[j][j]**0.5 )
    else:
        errs = [ np.inf, np.inf ]
    return np.append( pars, np.array( errs ) )

regr1 = linearRegression( xs, ys )
regr2 = linearRegressionSP( xs, ys )

print( regr1 )
print( 'Calculated sigma = %f' %  ( regr1[3] * np.sqrt( xs.shape[0] ) ) )
print( regr2 )
#print( 'B = %f must be in ( %f,\t%f )' % ( B, regr1[1] - regr1[3], regr1[1] + regr1[3] ) )


plt.plot( xs, ys, 'bo' )
plt.plot( xs, regr1[0] * xs + regr1[1] )
plt.show()

输出结果：

[1.3481681543925064, 2.1729338701374137, 0.0036028493647274687, 0.0062446292528624348]
Calculated sigma = 0.237624 # quite close to yerr
[ 1.34816815  2.17293387  0.00360534  0.01907908]

很有意思的是，curvefit和bootstrap会给出什么样的结果呢……

Answer 2

我在寻找自己类似问题的答案时发现了你的问题。简单来说，leastsq 输出的 cov_x 需要乘以残差方差，也就是：

s_sq = (func(popt, args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq

就像在 curve_fit.py 中那样。这是因为 leastsq 输出的是一个分数协方差矩阵。我遇到的一个大问题是，残差方差在网上搜索时显示的内容不太一样。

其实，残差方差就是你拟合结果中的简化卡方值。

Answer 3

更新于2016年4月6日

在拟合参数中获取正确的误差可能会很微妙。

我们来想象一下拟合一个函数 y=f(x)，你有一组数据点 (x_i, y_i, yerr_i)，其中 i 是每个数据点的索引。

在大多数物理测量中，误差 yerr_i 通常是测量设备或过程的系统性不确定性，因此可以认为它是一个常数，不依赖于 i。

该使用哪个拟合函数，以及如何获得参数误差？

optimize.leastsq 方法会返回一个分数协方差矩阵。将这个矩阵的所有元素乘以残差方差（即减少的卡方值），然后对对角线元素开平方，就可以得到拟合参数的标准差估计。我在下面的某个函数中包含了实现这个的代码。

另一方面，如果你使用 optimize.curvefit，上述过程的第一部分（乘以减少的卡方值）会在后台为你完成。然后你只需要对协方差矩阵的对角线元素开平方，就能得到拟合参数的标准差估计。

此外，optimize.curvefit 提供了可选参数，以处理更一般的情况，其中每个数据点的 yerr_i 值是不同的。有关更多信息，可以查看 文档：

sigma : None or M-length sequence, optional
    If not None, the uncertainties in the ydata array. These are used as
    weights in the least-squares problem
    i.e. minimising ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )``
    If None, the uncertainties are assumed to be 1.

absolute_sigma : bool, optional
    If False, `sigma` denotes relative weights of the data points.
    The returned covariance matrix `pcov` is based on *estimated*
    errors in the data, and is not affected by the overall
    magnitude of the values in `sigma`. Only the relative
    magnitudes of the `sigma` values matter.

我怎么能确保我的误差是正确的？

确定拟合参数的标准误差的正确估计是一个复杂的统计问题。optimize.curvefit 和 optimize.leastsq 实现的协方差矩阵结果实际上依赖于对误差的概率分布和参数之间相互作用的假设；这些相互作用可能存在，具体取决于你的拟合函数 f(x)。

在我看来，处理复杂的 f(x) 最好的方法是使用自助法，这在 这个链接 中有详细说明。

让我们看看一些例子

首先，一些基础代码。我们定义一个波动的函数，并生成一些带有随机误差的数据。我们将生成一个具有小随机误差的数据集。

import numpy as np
from scipy import optimize
import random

def f( x, p0, p1, p2):
    return p0*x + 0.4*np.sin(p1*x) + p2

def ff(x, p):
    return f(x, *p)

# These are the true parameters
p0 = 1.0
p1 = 40
p2 = 2.0

# These are initial guesses for fits:
pstart = [
    p0 + random.random(),
    p1 + 5.*random.random(), 
    p2 + random.random()
]

%matplotlib inline
import matplotlib.pyplot as plt
xvals = np.linspace(0., 1, 120)
yvals = f(xvals, p0, p1, p2)

# Generate data with a bit of randomness
# (the noise-less function that underlies the data is shown as a blue line)

xdata = np.array(xvals)
np.random.seed(42)
err_stdev = 0.2
yvals_err =  np.random.normal(0., err_stdev, len(xdata))
ydata = f(xdata, p0, p1, p2) + yvals_err

plt.plot(xvals, yvals)
plt.plot(xdata, ydata, 'o', mfc='None')

现在，让我们使用各种可用的方法来拟合这个函数：

`optimize.leastsq`

def fit_leastsq(p0, datax, datay, function):

    errfunc = lambda p, x, y: function(x,p) - y

    pfit, pcov, infodict, errmsg, success = \
        optimize.leastsq(errfunc, p0, args=(datax, datay), \
                          full_output=1, epsfcn=0.0001)

    if (len(datay) > len(p0)) and pcov is not None:
        s_sq = (errfunc(pfit, datax, datay)**2).sum()/(len(datay)-len(p0))
        pcov = pcov * s_sq
    else:
        pcov = np.inf

    error = [] 
    for i in range(len(pfit)):
        try:
          error.append(np.absolute(pcov[i][i])**0.5)
        except:
          error.append( 0.00 )
    pfit_leastsq = pfit
    perr_leastsq = np.array(error) 
    return pfit_leastsq, perr_leastsq 

pfit, perr = fit_leastsq(pstart, xdata, ydata, ff)

print("\n# Fit parameters and parameter errors from lestsq method :")
print("pfit = ", pfit)
print("perr = ", perr)

# Fit parameters and parameter errors from lestsq method :
pfit =  [  1.04951642  39.98832634   1.95947613]
perr =  [ 0.0584024   0.10597135  0.03376631]

`optimize.curve_fit`

def fit_curvefit(p0, datax, datay, function, yerr=err_stdev, **kwargs):
    """
    Note: As per the current documentation (Scipy V1.1.0), sigma (yerr) must be:
        None or M-length sequence or MxM array, optional
    Therefore, replace:
        err_stdev = 0.2
    With:
        err_stdev = [0.2 for item in xdata]
    Or similar, to create an M-length sequence for this example.
    """
    pfit, pcov = \
         optimize.curve_fit(f,datax,datay,p0=p0,\
                            sigma=yerr, epsfcn=0.0001, **kwargs)
    error = [] 
    for i in range(len(pfit)):
        try:
          error.append(np.absolute(pcov[i][i])**0.5)
        except:
          error.append( 0.00 )
    pfit_curvefit = pfit
    perr_curvefit = np.array(error)
    return pfit_curvefit, perr_curvefit 

pfit, perr = fit_curvefit(pstart, xdata, ydata, ff)

print("\n# Fit parameters and parameter errors from curve_fit method :")
print("pfit = ", pfit)
print("perr = ", perr)

# Fit parameters and parameter errors from curve_fit method :
pfit =  [  1.04951642  39.98832634   1.95947613]
perr =  [ 0.0584024   0.10597135  0.03376631]

`bootstrap`

def fit_bootstrap(p0, datax, datay, function, yerr_systematic=0.0):

    errfunc = lambda p, x, y: function(x,p) - y

    # Fit first time
    pfit, perr = optimize.leastsq(errfunc, p0, args=(datax, datay), full_output=0)


    # Get the stdev of the residuals
    residuals = errfunc(pfit, datax, datay)
    sigma_res = np.std(residuals)

    sigma_err_total = np.sqrt(sigma_res**2 + yerr_systematic**2)

    # 100 random data sets are generated and fitted
    ps = []
    for i in range(100):

        randomDelta = np.random.normal(0., sigma_err_total, len(datay))
        randomdataY = datay + randomDelta

        randomfit, randomcov = \
            optimize.leastsq(errfunc, p0, args=(datax, randomdataY),\
                             full_output=0)

        ps.append(randomfit) 

    ps = np.array(ps)
    mean_pfit = np.mean(ps,0)

    # You can choose the confidence interval that you want for your
    # parameter estimates: 
    Nsigma = 1. # 1sigma gets approximately the same as methods above
                # 1sigma corresponds to 68.3% confidence interval
                # 2sigma corresponds to 95.44% confidence interval
    err_pfit = Nsigma * np.std(ps,0) 

    pfit_bootstrap = mean_pfit
    perr_bootstrap = err_pfit
    return pfit_bootstrap, perr_bootstrap 

pfit, perr = fit_bootstrap(pstart, xdata, ydata, ff)

print("\n# Fit parameters and parameter errors from bootstrap method :")
print("pfit = ", pfit)
print("perr = ", perr)

# Fit parameters and parameter errors from bootstrap method :
pfit =  [  1.05058465  39.96530055   1.96074046]
perr =  [ 0.06462981  0.1118803   0.03544364]

观察

我们已经开始看到一些有趣的事情，所有三种方法的参数和误差估计几乎一致。这很好！

现在，假设我们想告诉拟合函数，我们的数据中还有其他的不确定性，可能是一个系统性的不确定性，会导致额外的误差，值是 err_stdev 的二十倍。这是 很多 的误差，实际上，如果我们模拟一些具有这种误差的数据，它会看起来像这样：

在这种噪声水平下，我们显然没有希望恢复拟合参数。

首先，让我们意识到 leastsq 甚至不允许我们输入这个新的系统性误差信息。让我们看看 curve_fit 在我们告诉它误差时会发生什么：

pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev)

print("\nFit parameters and parameter errors from curve_fit method (20x error) :")
print("pfit = ", pfit)
print("perr = ", perr)

Fit parameters and parameter errors from curve_fit method (20x error) :
pfit =  [  1.04951642  39.98832633   1.95947613]
perr =  [ 0.0584024   0.10597135  0.03376631]

什么？？这肯定是错的！

这本来是故事的结尾，但最近 curve_fit 添加了 absolute_sigma 可选参数：

pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev, absolute_sigma=True)

print("\n# Fit parameters and parameter errors from curve_fit method (20x error, absolute_sigma) :")
print("pfit = ", pfit)
print("perr = ", perr)

# Fit parameters and parameter errors from curve_fit method (20x error, absolute_sigma) :
pfit =  [  1.04951642  39.98832633   1.95947613]
perr =  [ 1.25570187  2.27847504  0.72600466]

这稍微好一些，但仍然有点可疑。curve_fit 认为我们可以从这个嘈杂的信号中得到一个拟合，p1 参数的误差水平为10%。让我们看看 bootstrap 的结果：

pfit, perr = fit_bootstrap(pstart, xdata, ydata, ff, yerr_systematic=20.0)

print("\nFit parameters and parameter errors from bootstrap method (20x error):")
print("pfit = ", pfit)
print("perr = ", perr)

Fit parameters and parameter errors from bootstrap method (20x error):
pfit =  [  2.54029171e-02   3.84313695e+01   2.55729825e+00]
perr =  [  6.41602813  13.22283345   3.6629705 ]

啊，这可能是对我们拟合参数误差的更好估计。bootstrap 认为它对 p1 的不确定性大约为34%。

总结

optimize.leastsq 和 optimize.curvefit 为我们提供了一种估计拟合参数误差的方法，但我们不能在使用这些方法时不加思考。bootstrap 是一种使用暴力的方法，在我看来，它在一些更难以解释的情况下往往效果更好。

我强烈建议你查看一个具体的问题，并尝试 curvefit 和 bootstrap。如果它们相似，那么 curvefit 计算起来更便宜，值得使用。如果它们差异显著，那么我会选择 bootstrap。