<p>2016年4月6日更新</p>
<h2>在大多数情况下,在拟合参数中获得正确的错误是很微妙的。</h2>
<p>让我们考虑拟合一个函数<code>y=f(x)</code>,对于这个函数,您有一组数据点<code>(x_i, y_i, yerr_i)</code>,其中<code>i</code>是在每个数据点上运行的索引。</p>
<p>在大多数物理测量中,误差<code>yerr_i</code>是测量装置或程序的系统不确定度,因此可以认为它是一个不依赖于<code>i</code>的常数。</p>
<h3>使用哪种拟合函数,如何获得参数误差?</h3>
<p><code>optimize.leastsq</code>方法将返回分数协方差矩阵。将该矩阵的所有元素乘以剩余方差(即减小的卡方)并取对角元素的平方根,将得到拟合参数的标准偏差估计值。我已经在下面的一个函数中包含了这样做的代码。</p>
<p>另一方面,如果使用<code>optimize.curvefit</code>,则上述过程的第一部分(乘以缩减的卡方)是在幕后为您完成的。然后需要取协方差矩阵对角线元素的平方根,以获得拟合参数的标准偏差估计值。</p>
<p>此外,<code>optimize.curvefit</code>提供可选参数来处理更一般的情况,其中每个数据点的<code>yerr_i</code>值不同。从<a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html" rel="noreferrer">documentation</a>:</p>
<pre><code>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.
</code></pre>
<h3>我怎样才能确定我的错误是正确的?</h3>
<p>确定合适的拟合参数标准误差估计是一个复杂的统计问题。由<code>optimize.curvefit</code>和<code>optimize.leastsq</code>实现的协方差矩阵的结果实际上依赖于关于误差概率分布和参数之间相互作用的假设;可能存在的相互作用,取决于特定的拟合函数<code>f(x)</code>。</p>
<p>在我看来,处理复杂的<code>f(x)</code>的最佳方法是使用bootstrap方法,如<a href="http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html" rel="noreferrer">this link</a>所述。</p>
<h3>让我们看一些例子</h3>
<p>首先,一些样板代码。让我们定义一个曲线函数并生成一些带有随机错误的数据。我们将生成一个带有小随机误差的数据集。</p>
<pre><code>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')
</code></pre>
<p><a href="https://i.stack.imgur.com/msQU1.png" rel="noreferrer"><img src="https://i.stack.imgur.com/msQU1.png" alt="fig01"/></a></p>
<p>现在,让我们使用各种可用的方法来拟合函数:</p>
<h3>`优化.租赁
<pre><code>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)
</code></pre>
<p><br/></p>
<pre><code># Fit parameters and parameter errors from lestsq method :
pfit = [ 1.04951642 39.98832634 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
</code></pre>
<p><br/></p>
<h3>`优化.曲线拟合`</h3>
<pre><code>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)
</code></pre>
<p><br/></p>
<pre><code># Fit parameters and parameter errors from curve_fit method :
pfit = [ 1.04951642 39.98832634 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
</code></pre>
<p><br/></p>
<h3>`引导`</h3>
<pre><code>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)
</code></pre>
<p><br/></p>
<pre><code># Fit parameters and parameter errors from bootstrap method :
pfit = [ 1.05058465 39.96530055 1.96074046]
perr = [ 0.06462981 0.1118803 0.03544364]
</code></pre>
<p><br/></p>
<h2>观察</h2>
<p>我们已经开始看到一些有趣的东西,三种方法的参数和误差估计几乎一致。太好了!</p>
<p>现在,假设我们想告诉拟合函数,在我们的数据中还有一些其他的不确定性,可能是一个系统性的不确定性,它会造成20倍于<code>err_stdev</code>值的额外误差。实际上,如果我们用这种错误来模拟一些数据,它看起来会是这样的:</p>
<p><a href="https://i.stack.imgur.com/RRljw.png" rel="noreferrer"><img src="https://i.stack.imgur.com/RRljw.png" alt="enter image description here"/></a></p>
<p>我们当然不希望用这种噪声水平来恢复拟合参数。</p>
<p>首先,让我们认识到<code>leastsq</code>甚至不允许我们输入这个新的系统错误信息。让我们看看<code>curve_fit</code>当我们告诉它错误时会做什么:</p>
<pre><code>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)
</code></pre>
<p><br/></p>
<pre><code>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]
</code></pre>
<p>什么??这肯定是错的!</p>
<p>这曾经是故事的结尾,但最近<code>curve_fit</code>添加了<code>absolute_sigma</code>可选参数:</p>
<pre><code>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)
</code></pre>
<p><br/></p>
<pre><code># 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]
</code></pre>
<p>这有点好,但还是有点可疑。<code>curve_fit</code>认为我们可以在<code>p1</code>参数中以10%的误差拟合出噪声信号。让我们看看<code>bootstrap</code>要说什么:</p>
<pre><code>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)
</code></pre>
<p><br/></p>
<pre><code>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 ]
</code></pre>
<p>啊,这也许是对拟合参数误差的更好估计。<code>bootstrap</code>认为它知道<code>p1</code>,不确定度约为34%。</p>
<h2>摘要</h2>
<p><code>optimize.leastsq</code>和<code>optimize.curvefit</code>为我们提供了一种估计拟合参数误差的方法,但我们不能只使用这些方法而不稍微质疑它们。<code>bootstrap</code>是一种使用暴力的统计方法,在我看来,它有一种在可能更难解释的情况下更好地工作的趋势。</p>
<p>我强烈建议查看一个特定的问题,并尝试<code>curvefit</code>和<code>bootstrap</code>。如果它们相似,那么<code>curvefit</code>计算起来要便宜得多,因此可能值得使用。如果它们有显著的不同,那么我的钱就在<code>bootstrap</code>上。</p>