使用 python 中的 optimize.leastsq 方法获取拟合参数的标准错误
Posted
技术标签:
【中文标题】使用 python 中的 optimize.leastsq 方法获取拟合参数的标准错误【英文标题】：Getting standard errors on fitted parameters using the optimize.leastsq method in python 【发布时间】：2013-01-12 22:22:14 【问题描述】：我有一组数据（位移与时间），我使用 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 是时间和位移值的数组（基本上只有 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]]
我怀疑目前的合身性有点可疑。这将在我可以解决错误时得到确认。
【问题讨论】：
docs.scipy.org/doc/scipy/reference/generated/… 不完全确定您的要求，也许一些例子会有所帮助。 这是我一直在看的主要文档。我使用的代码已经添加到上面的描述中 我想我已经找到了一种解决方法（尽管在重写代码方面有点不方便）我认为“optimise.curve_fit”输出协方差矩阵，您可以从中获取错误，它使用与“optimize.leastsq”相同的最小二乘回归方法。有人可以确认这是正确的吗？ 是的，curve_fit 返回参数估计的协方差矩阵（不确定性）。如果想直接使用leastsq，也可以查看curve_fit的来源。 【参考方案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 值不同。来自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.
如何确定我的错误是正确的？
确定拟合参数中标准误差的正确估计是一个复杂的统计问题。由optimize.curvefit 和optimize.leastsq 实现的协方差矩阵的结果实际上依赖于关于错误概率分布和参数之间相互作用的假设；可能存在的交互，具体取决于您的特定拟合函数f(x)。
在我看来，处理复杂的f(x) 的最佳方法是使用引导方法，该方法在this link 中进行了概述。
让我们看一些例子
首先，一些样板代码。让我们定义一个波浪线函数并生成一些带有随机误差的数据。我们将生成一个随机误差较小的数据集。
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]
`引导程序`
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 值的 20 倍。这是很多的错误，事实上，如果我们模拟一些具有这种错误的数据，它会看起来像这样：
在这种噪音水平下，我们当然不可能恢复拟合参数。
首先，让我们意识到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 上。
【讨论】：
那么，让我看看我是否明白你做了什么。您根据提供的错误（或残差）围绕提供的Y 生成随机值，您执行此操作 100 次，在每个步骤中获取拟合参数。然后您将STD视为错误？与statsmodels的WLS相比如何？
嗨 Yotam，我不熟悉 statsmodels 中的 WLS 函数。我的目标是处理单独具有较大误差/不确定性但其平均值相当适合模型的测量。目标是提供拟合参数误差的合理值。即使平均值完全位于拟合线上，如果每次测量的误差都很大，那么拟合参数误差也必须反映这一点。
一种方法是“蛮力”引导方法（我在这里展示过），另一种方法是研究误差函数的依赖性（卡方，$\chi^2 $) 在每个拟合参数上，并为每个参数建立一个置信区间。有时人们使用误差函数对拟合参数 ($\mathrmd\chi^2 / \mathrmdp$ ) 的已知导数来快速访问这个置信区间，或者他们在数值上改变所有其他参数保持固定的参数，并根据观察到的 $\chi^2$ 的变化定义置信区间
@ThePredator 默认情况下，sigma 值应仅被理解为单个数据点的权重。 curvefit 将使用权重，并将对拟合值产生影响，但对协方差矩阵不。当您说 absolute_sigma=True 时，您是在告诉 curvefit 您的 sigma 不仅是权重，而且是与 $y$ 相同单位的实际不确定性。在这种情况下，curvefit 继续并在pcov 中返回未归一化的协方差矩阵。这样就可以对pcov中的对角线项取平方根，直接得到参数误差。
@PalSzabo 我刚刚编辑了一个文档字符串来解释如何在示例中解决该问题 :)【参考方案2】：
在线性回归的情况下可以精确计算误差。实际上，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 给出的结果很有趣...
【讨论】：
【参考方案3】：在尝试回答我自己的类似question 时发现了您的问题。
简短的回答。最小平方输出的cov_x 应该乘以残差方差。即
s_sq = (func(popt, args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq
如curve_fit.py。这是因为 minimumsq 输出分数协方差矩阵。我最大的问题是，在谷歌搜索时，剩余方差会显示为其他内容。
残差只是从您的拟合中减少卡方。
【讨论】：
以上是关于使用 python 中的 optimize.leastsq 方法获取拟合参数的标准错误的主要内容，如果未能解决你的问题，请参考以下文章