python optimize_使用python中的optimize.leastsq方法获取拟合参数的标准错误

于4/6/2016更新

在大多数情况下，在拟合参数中获得正确的错误可能是微妙的。

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

在大多数物理测量中，误差yerr_i是测量装置或程序的系统不确定性，因此可以认为是不依赖于i的常数。

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

optimize.leastsq方法将返回分数协方差矩阵。将该矩阵的所有元素乘以剩余方差(即减小的奇偶平方)并取乘角元素的平方根将给出拟合参数的标准偏差的估计。我已经将代码包含在下面的一个函数中。

另一方面，如果您使用optimize.curvefit，则会在程序中为您执行上述步骤的第一部分(乘以减少的卡方)。在这种情况下，您只需要取协方差矩阵的对角元素的平方根来得到拟合参数的标准偏差的估计值。

此外，optimize.curvefit提供了可选参数来处理更一般的情况，其中yerr_i值不是独立于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

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

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("\nFit paramters 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):

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("\nFit 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("\nFit 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值的二十倍的额外的误差。这是一个很大的错误，事实上，如果我们模拟一些数据与这种错误，它将如下所示：

我们肯定没有希望，我们可以用这种噪音水平恢复适合的参数。

首先，我们认识到，最不发达国家甚至不允许我们输入这个新的系统错误信息。当我们告诉它有关错误的时候，我们来看看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("\nFit 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为我们提供了一种估计拟合参数中的错误的方法，但是我们不能仅仅使用这些方法而不用质疑它们。引导是一种使用暴力的统计方法，在我看来，它有可能在更难解释的情况下更好地工作。

我强烈推荐看一个特定的问题，并尝试曲线和引导。如果它们相似，那么曲线拟合计算便宜得多，所以很可能使用。如果他们有显着差异，那么我的钱就会在引导下。