使用 python 中的 optimize.leastsq 方法获取拟合参数的标准误差 [英] Getting standard errors on fitted parameters using the optimize.leastsq method in python

问题描述

我有一组数据(位移 vs 时间)，我使用 optimize.leastsq 方法将这些数据拟合到几个方程中.我现在正在寻找拟合参数的错误值.查看文档，输出的矩阵是雅可比矩阵，我必须将其乘以残差矩阵才能得到我的值.不幸的是，我不是统计学家，所以我对术语有些不知所措.

据我所知，我需要的是与我的拟合参数相匹配的协方差矩阵，因此我可以对对角元素进行平方根以得到拟合参数的标准误差.我有一个模糊的阅读记忆，协方差矩阵无论如何都是从 optimize.leastsq 方法输出的.这样对吗?如果不是，你将如何让残差矩阵与输出的雅可比矩阵相乘以获得我的协方差矩阵?

任何帮助将不胜感激.我对 python 很陌生，因此如果问题是一个基本问题，我深表歉意.

拟合代码如下:

fitfunc = lambda p, t: p[0]+p[1]*np.log(t-p[2])+ p[3]*t #目标函数'errfunc = lambda p, t, y: (fitfunc(p, t) - y)# 到目标函数的距离p0 = [ 1,1,1,1] # 参数的初始猜测out = optimize.leastsq(errfunc, p0[:], args=(t, disp,), full_output=1)

args t 和 disp 是时间和位移值数组(基本上只有 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]]

无论如何，我怀疑现在的合身有点可疑.当我可以解决错误时，这将得到确认.

解决方案

更新于 4/6/2016

在大多数情况下，获取拟合参数中的正确错误可能很微妙.

让我们考虑拟合一个函数 y=f(x)，其中有一组数据点 (x_i, y_i, yerr_i)，其中 i 是运行在每个数据点上的索引.

在大多数物理测量中，误差 yerr_i 是测量设备或程序的系统不确定性，因此可以将其视为不依赖于 i.

使用哪个拟合函数，以及如何获取参数误差?

optimize.leastsq 方法将返回分数协方差矩阵.将此矩阵的所有元素乘以残差方差(即减少的卡方)并取对角元素的平方根，您将得到拟合参数标准偏差的估计值.我已经在以下功能之一中包含了执行此操作的代码.

另一方面，如果您使用optimize.curvefit，则上述过程的第一部分(乘以减少的卡方)是在幕后为您完成的.然后，您需要取协方差矩阵的对角元素的平方根，以得到拟合参数标准差的估计值.

此外，optimize.curvefit 提供可选参数来处理更一般的情况，其中每个数据点的 yerr_i 值不同.来自

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

`optimize.leastsq`

def fit_leastsq(p0, datax, datay, function):errfunc = lambda p, x, y: function(x,p) - ypfit, pcov, infodict, errmsg, 成功 = \optimize.leastsq(errfunc, p0, args=(datax, datay), \full_output=1, epsfcn=0.0001)如果 (len(datay) > len(p0)) 并且 pcov 不是 None:s_sq = (errfunc(pfit, datax, datay)**2).sum()/(len(datay)-len(p0))pcov = pcov * s_sq别的:pcov = np.inf错误 = []对于范围内的 i(len(pfit)):尝试:error.append(np.absolute(pcov[i][i])**0.5)除了:错误.追加(0.00)pfit_leastsq = pfitperr_leastsq = np.array(error)返回 pfit_leastsq、perr_leastsqpfit, perr = fit_leastsq(pstart, xdata, ydata, ff)print("\n#拟合参数和来自lestsq方法的参数错误:")打印(pfit ="，pfit)打印("perr = ", perr)




# 从 lestsq 方法中拟合参数和参数错误: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):"""注意:根据当前文档(Scipy V1.1.0)，sigma (yerr) 必须是:无或 M 长度序列或 MxM 数组，可选因此，替换:err_stdev = 0.2和:err_stdev = [xdata 中的项目为 0.2]或者类似的，为这个例子创建一个 M 长度的序列."""pfit，pcov = \optimize.curve_fit(f,datax,datay,p0=p0,\sigma=yerr, epsfcn=0.0001, **kwargs)错误 = []对于范围内的 i(len(pfit)):尝试:error.append(np.absolute(pcov[i][i])**0.5)除了:错误.追加(0.00)pfit_curvefit = pfitperr_curvefit = np.array(error)返回 pfit_curvefit, perr_curvefitpfit, perr = fit_curvefit(pstart, xdata, ydata, ff)print("\n#从curve_fit方法中拟合参数和参数错误:")打印(pfit ="，pfit)打印("perr = ", perr)




# 拟合参数和来自curve_fit方法的参数错误:pfit = [ 1.04951642 39.98832634 1.95947613]perr = [ 0.0584024 0.10597135 0.03376631]




`引导程序`

def fit_bootstrap(p0, datax, datay, function, yerr_systematic=0.0):errfunc = lambda p, x, y: function(x,p) - y# 第一次拟合pfit, perr = optimize.leastsq(errfunc, p0, args=(datax, datay), full_output=0)# 获取残差的标准差残差 = errfunc(pfit, datax, datay)sigma_res = np.std(残差)sigma_err_total = np.sqrt(sigma_res**2 + yerr_systematic**2)# 生成并拟合 100 个随机数据集ps = []对于我在范围内(100):randomDelta = np.random.normal(0., sigma_err_total, len(datay))randomdataY = datay + randomDelta随机拟合，randomcov = \optimize.leastsq(errfunc, p0, args=(datax, randomdataY),\full_output=0)ps.append(randomfit)ps = np.array(ps)mean_pfit = np.mean(ps,0)# 你可以选择你想要的置信区间# 参数估计:Nsigma = 1. # 1sigma 与上述方法大致相同# 1sigma 对应 68.3% 置信区间# 2sigma 对应 95.44% 置信区间err_pfit = Nsigma * np.std(ps,0)pfit_bootstrap = mean_pfitperr_bootstrap = err_pfit返回 pfit_bootstrap、perr_bootstrappfit, perr = fit_bootstrap(pstart, xdata, ydata, ff)print("\n#从引导方法中拟合参数和参数错误:")打印(pfit ="，pfit)打印("perr = ", perr)




# 从引导方法中拟合参数和参数错误:pfit = [ 1.05058465 39.96530055 1.96074046]perr = [ 0.06462981 0.1118803 0.03544364]




观察

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

现在，假设我们要告诉拟合函数我们的数据中存在其他一些不确定性，可能是系统性不确定性，它会导致 err_stdev 值的 20 倍的额外误差.这是很多的错误，事实上，如果我们用这种错误模拟一些数据，它看起来像这样:

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

首先，让我们意识到leastsq 甚至不允许我们输入这个新的系统错误信息.让我们看看当我们告诉它错误时 curve_fit 做了什么:

pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev)print("\n曲线拟合方法中的拟合参数和参数错误(20x 错误):")打印(pfit ="，pfit)打印("perr = ", perr)




Fit 参数和来自 curve_fit 方法的参数错误(20x 错误):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# 拟合参数和曲线拟合方法的参数错误(20x 错误，absolute_sigma):")打印(pfit ="，pfit)打印("perr = ", perr)




# 拟合参数和来自curve_fit方法的参数误差(20x误差，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("\n从引导方法中拟合参数和参数错误(20x 错误):")打印(pfit ="，pfit)打印("perr = ", perr)




从引导方法中拟合参数和参数错误(20x 错误):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.

I have a set of data (displacement vs time) which I have fitted to a couple of equations using the optimize.leastsq method. I am now looking to get error values on the fitted parameters. Looking through the documentation the matrix outputted is the jacobian matrix, and I must multiply this by the residual matrix to get my values. Unfortunately I am not a statistician so I am drowning somewhat in the terminology.

From what I understand all I need is the covariance matrix that goes with my fitted parameters, so I can square root the diagonal elements to get my standard error on the fitted parameters. I have a vague memory of reading that the covariance matrix is what is output from the optimize.leastsq method anyway. Is this correct? If not how would you go about getting the residual matrix to multiply the outputted Jacobian by to get my covariance matrix?

Any help would be greatly appreciated. I am very new to python and therefore apologise if the question turns out to be a basic one.

the fitting code is as follows:

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)

The args t and disp is and array of time and displcement values (basically just 2 columns of data). I have imported everything needed at the tope of the code. The fitted values and the matrix provided by the output is as follows:

[  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]]

I suspect the fit is a little suspect anyway at the moment. This will be confirmed when I can get the errors out.

解决方案

Updated on 4/6/2016

Getting the correct errors in the fit parameters can be subtle in most cases.

Let's think about fitting a function y=f(x) for which you have a set of data points (x_i, y_i, yerr_i), where i is an index that runs over each of your data points.

In most physical measurements, the error yerr_i is a systematic uncertainty of the measuring device or procedure, and so it can be thought of as a constant that does not depend on i.

Which fitting function to use, and how to obtain the parameter errors?

The optimize.leastsq method will return the fractional covariance matrix. Multiplying all elements of this matrix by the residual variance (i.e. the reduced chi squared) and taking the square root of the diagonal elements will give you an estimate of the standard deviation of the fit parameters. I have included the code to do that in one of the functions below.

On the other hand, if you use optimize.curvefit, the first part of the above procedure (multiplying by the reduced chi squared) is done for you behind the scenes. You then need to take the square root of the diagonal elements of the covariance matrix to get an estimate of the standard deviation of the fit parameters.

Furthermore, optimize.curvefit provides optional parameters to deal with more general cases, where the yerr_i value is different for each data point. From the documentation:

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.

How can I be sure that my errors are correct?

Determining a proper estimate of the standard error in the fitted parameters is a complicated statistical problem. The results of the covariance matrix, as implemented by optimize.curvefit and optimize.leastsq actually rely on assumptions regarding the probability distribution of the errors and the interactions between parameters; interactions which may exist, depending on your specific fit function f(x).

In my opinion, the best way to deal with a complicated f(x) is to use the bootstrap method, which is outlined in this link.

Let's see some examples

First, some boilerplate code. Let's define a squiggly line function and generate some data with random errors. We will generate a dataset with a small random error.

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')

Now, let's fit the function using the various methods available:

`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]

Observations

We already start to see something interesting, the parameters and error estimates for all three methods nearly agree. That is good!

Now, suppose we want to tell the fitting functions that there is some other uncertainty in our data, perhaps a systematic uncertainty that would contribute an additional error of twenty times the value of err_stdev. That is a lot of error, in fact, if we simulate some data with that kind of error it would look like this:

There is certainly no hope that we could recover the fit parameters with this level of noise.

To begin with, let's realize that leastsq does not even allow us to input this new systematic error information. Let's see what curve_fit does when we tell it about the error:

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?? This must certainly be wrong!

This used to be the end of the story, but recently curve_fit added the absolute_sigma optional parameter:

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]

That is somewhat better, but still a little fishy. curve_fit thinks we can get a fit out of that noisy signal, with a level of 10% error in the p1 parameter. Let's see what bootstrap has to say:

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 ]

Ah, that is perhaps a better estimate of the error in our fit parameter. bootstrap thinks it knows p1 with about a 34% uncertainty.

Summary

optimize.leastsq and optimize.curvefit provide us a way to estimate errors in fitted parameters, but we cannot just use these methods without questioning them a little bit. The bootstrap is a statistical method which uses brute force, and in my opinion, it has a tendency of working better in situations that may be harder to interpret.

I highly recommend looking at a particular problem, and trying curvefit and bootstrap. If they are similar, then curvefit is much cheaper to compute, so probably worth using. If they differ significantly, then my money would be on the bootstrap.

这篇关于使用 python 中的 optimize.leastsq 方法获取拟合参数的标准误差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！