使用 python 中的 optimize.leastsq 方法获取拟合参数的标准误差
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/14581358/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Getting standard errors on fitted parameters using the optimize.leastsq method in python
提问by Phil
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.
我有一组数据（位移 vs 时间），我使用 optimize.leastsq 方法将它们拟合到几个方程中。我现在正在寻找拟合参数的错误值。查看文档，输出的矩阵是雅可比矩阵，我必须将其乘以残差矩阵才能得到我的值。不幸的是，我不是统计学家，所以我对术语有些不知所措。
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?
据我所知，我需要的是与我的拟合参数相匹配的协方差矩阵，因此我可以对对角线元素求平方以获得拟合参数的标准误差。我有一个模糊的阅读记忆，协方差矩阵无论如何都是从 optimize.leastsq 方法输出的。这样对吗？如果不是，你将如何让残差矩阵与输出的雅可比矩阵相乘以获得我的协方差矩阵？
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.
任何帮助将不胜感激。我对 python 很陌生，因此如果问题是一个基本问题，我深表歉意。
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:
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]]
I suspect the fit is a little suspect anyway at the moment. This will be confirmed when I can get the errors out.
无论如何，我怀疑现在的合身有点可疑。当我可以解决错误时，这将得到确认。
回答by HansHarhoff
Found your question while trying to answer my own similar question.
Short answer. The cov_xthat leastsq outputs should be multiplied by the residual variance. i.e.
在尝试回答我自己的类似问题时发现了您的问题。简短的回答。在cov_x该leastsq输出应当由剩余方差相乘。IE
s_sq = (func(popt, args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq
as in curve_fit.py. This is because leastsq outputs the fractional covariance matrix. My big problem was that residual variance shows up as something else when googling it.
如curve_fit.py. 这是因为 leastsq 输出分数协方差矩阵。我的大问题是残差方差在谷歌搜索时显示为其他东西。
Residual variance is simply reduced chi square from your fit.
残差只是从您的拟合中减少卡方。
回答by Pedro M Duarte
Updated on 4/6/2016
2016 年 4 月 6 日更新
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 iis an index that runs over each of your data points.
让我们考虑拟合一个y=f(x)具有一组数据点的函数(x_i, y_i, yerr_i)，其中i是一个在每个数据点上运行的索引。
In most physical measurements, the error yerr_iis 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.
在大多数物理测量中，误差yerr_i是测量设备或程序的系统不确定性，因此可以将其视为不依赖于 的常数i。
Which fitting function to use, and how to obtain the parameter errors?
使用哪个拟合函数，以及如何获取参数误差？
The optimize.leastsqmethod 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.
该optimize.leastsq方法将返回分数协方差矩阵。将此矩阵的所有元素乘以残差方差（即减少的卡方）并取对角元素的平方根，您将得到拟合参数标准偏差的估计值。我已经在以下功能之一中包含了执行此操作的代码。
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.
另一方面，如果您使用optimize.curvefit，则上述过程的第一部分（乘以减少的卡方）是在幕后为您完成的。然后，您需要取协方差矩阵的对角元素的平方根，以获得拟合参数标准偏差的估计值。
Furthermore, optimize.curvefitprovides optional parameters to deal with more general cases, where the yerr_ivalue is different for each data point. From the documentation:
此外，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.
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.curvefitand optimize.leastsqactually 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).
确定拟合参数中标准误差的正确估计是一个复杂的统计问题。协方差矩阵的结果，由关于误差的概率分布和参数之间的相互作用的假设来实现，optimize.curvefit并且optimize.leastsq实际上依赖于这些假设；可能存在的相互作用，具体取决于您的特定拟合函数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.
在我看来，处理复杂问题的最佳方法f(x)是使用 bootstrap 方法，此链接中概述了该方法。
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`
`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`
`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 lotof error, in fact, if we simulate some data with that kind of error it would look like this:
现在，假设我们想告诉拟合函数，我们的数据中还有一些其他的不确定性，也许是一种系统性的不确定性，它会造成 20 倍的额外误差err_stdev。这是一个很大的错误，事实上，如果我们用这种错误模拟一些数据，它看起来像这样：
There is certainly no hope that we could recover the fit parameters with this level of noise.
我们当然没有希望在这种噪声水平下恢复拟合参数。
To begin with, let's realize that leastsqdoes not even allow us to input this new systematic error information. Let's see what curve_fitdoes when we tell it about the error:
首先，让我们意识到这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?? This must certainly be wrong!
什么？？这肯定是错误的！
This used to be the end of the story, but recently curve_fitadded the absolute_sigmaoptional parameter:
这曾经是故事的结尾，但最近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]
That is somewhat better, but still a little fishy. curve_fitthinks we can get a fit out of that noisy signal, with a level of 10% error in the p1parameter. Let's see what bootstraphas to say:
这样好一些，但仍然有点可疑。 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 ]
Ah, that is perhaps a better estimate of the error in our fit parameter. bootstrapthinks it knows p1with about a 34% uncertainty.
啊，这可能是对拟合参数误差的更好估计。bootstrap认为它知道p1大约 34% 的不确定性。
Summary
概括
optimize.leastsqand optimize.curvefitprovide us a way to estimate errors in fitted parameters, but we cannot just use these methods without questioning them a little bit. The bootstrapis 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.
optimize.leastsq并optimize.curvefit为我们提供了一种估计拟合参数误差的方法，但我们不能只使用这些方法而不会对它们提出一点质疑。这bootstrap是一种使用蛮力的统计方法，在我看来，它倾向于在可能难以解释的情况下更好地工作。
I highly recommend looking at a particular problem, and trying curvefitand bootstrap. If they are similar, then curvefitis much cheaper to compute, so probably worth using. If they differ significantly, then my money would be on the bootstrap.
我强烈建议查看一个特定的问题，然后尝试curvefit和bootstrap. 如果它们相似，则curvefit计算成本要低得多，因此可能值得使用。如果它们有显着差异，那么我的钱将放在bootstrap.
回答by jambulat
It is possible to calculate exactly the errors in the case of linear regression. And indeed the leastsq function gives the values which are different:
在线性回归的情况下，可以准确地计算误差。事实上，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()
output:
输出：
[1.3481681543925064, 2.1729338701374137, 0.0036028493647274687, 0.0062446292528624348]
Calculated sigma = 0.237624 # quite close to yerr
[ 1.34816815 2.17293387 0.00360534 0.01907908]
It's interesting which results curvefit and bootstrap will give...
有趣的是，curvefit 和 bootstrap 会给出哪些结果......
