使用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)

args t和disp是time和displcement值的数组(基本上只有2列数据).我已经导入了代码顶部所需的所有内容.输出提供的拟合值和矩阵如下:

[  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

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)是使用bootstrap方法,该方法在此链接中列出.

我们来看一些例子

首先,一些样板代码.让我们定义一个波浪线函数并生成一些随机错误的数据.我们将生成一个具有小随机误差的数据集.

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]

意见

我们已经开始看到一些有趣的东西,所有三种方法的参数和误差估计几乎都是一致的.那很好!

现在,假设我们想要告诉拟合函数我们的数据中存在一些其他不确定性,也许系统的不确定性将导致额外的误差为20倍的值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]

Whaat？这肯定是错的!

这曾经是故事的结尾,但最近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.

Answer 2

在尝试回答我自己的类似问题时发现了您的问题.简短的回答.在cov_x该leastsq输出应当由剩余方差相乘.即

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

如在curve_fit.py.这是因为leastsq输出分数协方差矩阵.我的一个大问题是,谷歌搜索时剩余差异显示为其他东西.

残余方差可以简化为适合您的卡方.