在数据分析过程中，我们经常会遇到线型拟合问题。标准的做法是用最小二乘法(least-squares method)来计算相关系数，用协方差矩阵(covariance matrix)估计误差。不过这是在假定所有数据点都是绝对精确的情况下才成立。而在实际工作中，测量值不可避免地带有误差，忽略这些误差显然会低估相关系数的误差。那么该如何合理考虑数据点的误差，并绘制相应的置信带(confidence band，表示拟合函数的可能出现范围)呢？这里我们借助Python中的相关函数来看一下。
上图中的数据座标和误差如下：
|
1 2 3 |
x = [1.2, 2.1, 2.9, 3.8, 4.9] # x轴位置 y = [1.1, 1.9, 3.0, 3.9, 5.1] # y轴位置 e = [0.3, 0.4, 0.5, 0.3, 0.4] # y轴误差 |
它们的相关关系可以用numpy模块中的多项式拟合函数 polyfit 来计算。
|
1 2 3 4 5 6 |
# fitting with polyfit import numpy as np popt, cov = np.polyfit(x,y,1, cov=True) # 1次多项式拟合，返回协方差矩阵 err = np.sqrt(np.diag(cov)) # 由协方差矩阵计算多项式系数误差 print(popt,err) # 输出拟合结果及误差 print(np.poly1d(popt)) # 输出拟合函数 |
可以得到拟合直线为 y = 1.10 x – 0.27， 由协方差矩阵得到的参数误差分别为 [0.03，0.11]。但是这个拟合并没有考虑数据点的误差。我们可以采用更灵活的拟合函数——scipy.optimize 模块中的非线性拟合函数 curve_fit 来实现，代码如下：
|
1 2 3 4 5 6 7 |
from scipy.optimize import curve_fit f = lambda x, *p: np.polyval(p, x) # 定义待拟合函数 # 依次传入拟合函数、拟合数据，参数初值，及误差 popt, cov = curve_fit(f, x, y, [1, 1], sigma=e, absolute_sigma = True) err = np.sqrt(np.diag(cov)) print(popt,err) # 输出拟合结果及误差 print(np.poly1d(popt)) # 输出拟合函数 |
可以得到拟合直线为 y = 1.09x-0.25，参数误差为[0.12，0.38]。拟合结果稍有不同，参数误差也合理增加。接下来就是根据这个拟合结果来绘制置信带。但置信带并非参数误差的简单叠加。我们可以测试一下，将误差直接与拟合结果叠加看效果如何，代码及结果如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 |
# manual check import pylab as plt xx = np.linspace(1,5,5) aa = popt[0] ; bb = popt[1] da = err[0] ; db = err[1] plt.errorbar(x,y,yerr=e, fmt="o",capsize=4) # 数据点 plt.plot(xx, np.poly1d([aa,bb])(xx),'r-') # 绘制最佳拟合结果 plt.plot(xx, np.poly1d([aa+da,bb+db])(xx),'c-') plt.plot(xx, np.poly1d([aa-da,bb-db])(xx),'c-') plt.plot(xx, np.poly1d([aa+da,bb-db])(xx),'m-') plt.plot(xx, np.poly1d([aa-da,bb+db])(xx),'m-') plt.show() |
可以看到这样得到的置信带并不正确。因为多项式系数之间并非完全独立(cov非对角阵)，它们的上下限自然无法同时达到。那么正确的做法是怎样的呢？我们可以重新看一下置信带的含义。它是表示真实函数以特定概率出现的区域。换句话说，1σ误差置信带就意味着真实函数有68%的可能性出现在这个区域。我们可以根据拟合结果模拟出足够多的函数，然后找出其中68%的结果集中的区域。实际操作时，我们统计的是模拟函数给出的模拟数据点的分布。具体代码和结果如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 |
xi = np.linspace(1, 5, 20) # 取样点越密，置信带越平滑 # 根据拟合结果生成1000个模拟函数 ps = np.random.multivariate_normal(popt, cov, 1000) # 根据模拟函数统计每个取样点上的y值 ysample = np.asarray([f(xi, *pi) for pi in ps]) # 在每个取样点上提取 68% 区间上下边界 lower = np.percentile(ysample, 16, axis=0) upper = np.percentile(ysample, 84, axis=0) plt.errorbar(x,y,yerr=e, fmt="o",capsize=4) plt.plot(xi, np.poly1d(popt)(xi), 'r-') # 绘制最佳拟合函数 plt.fill_between(xi, lower, upper, color="gray",alpha=0.2) plt.show() |
这个结果就很像样了。事实上，借助 uncertainties 模块我们可以由协方差矩阵直接得到置信带，代码如下：
|
1 2 3 4 5 6 7 8 9 |
import uncertainties as unc a, b = unc.correlated_values(popt, cov) # 将协方差矩阵信息导入系数 py = a * xi + b # 将采样点xi带入拟合函数得到目标数据特征 nom = unc.unumpy.nominal_values(py) # 得出各采样点处的期望值 std = unc.unumpy.std_devs(py) # 得出各采样点处的方差 plt.errorbar(x,y,yerr=e, fmt="o",capsize=4) # 绘制数据点 plt.plot(xi, nom, c='r') # 绘制最佳拟合函数 plt.plot(xi, nom - std, c='c') # 绘制最佳拟合的1σ下限 plt.plot(xi, nom + std, c='c') # 绘制最佳拟合的1σ上限 |
可以看到这样直接得到的青色置信带和前一个模拟方法得到的灰色置信带完全一致。对于通常的拟合需求来说，这样的结果已经可以用了。如果x轴位置上也有误差，那么curvefit就不再适用了，需要借助正交距离回归（Orthogonal Distance Regression，ODR）等方法。这里介绍一个更为普适通用的方法——自助法(bootstrap)。
既然每个数据点及其误差都代表一个高斯分布，我可以根据这些分布进行多次模拟采样，得到符合相同分布的多个模拟数据集。通过对这些模拟数据集的直接拟合，可以得到拟合函数的分布情况，进而求出其中位值和68.2%的置信带。具体代码和结果如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
func = lambda x, *p: np.polyval(p, x) n = 10000 # 生成模拟数据集 ysample = [] for i in range(len(y)): ysample.append(list(np.random.normal(y[i],e[i],n))) ysample = np.array(ysample) pset = [] yp = [] xi = np.linspace(1, 5, 20) for i in range(n): popt, cov = curve_fit(func, x, ysample[:,i], [1, 1])# 不带误差拟合 pset.append(popt) # 收集拟合参数，用于误差分析 yp = yp + [func(xi, *popt)] # 收集拟合结果点集，用于置信带绘制 # 确定参数及误差 pset = np.array(pset) par_median = np.percentile(pset, 50, axis=0) par_lower = np.percentile(pset, 15.9, axis=0) par_upper = np.percentile(pset, 84.1, axis=0) print("Best Fit:",par_median) print("Error:",par_median-par_lower, par_upper-par_median) # 确定置信带边界 yp = np.array(yp) l_fit = np.percentile(yp, 50, axis=0) l_lower = np.percentile(yp, 15.9, axis=0) l_upper = np.percentile(yp, 84.1, axis=0) l_lower2 = np.percentile(yp, 2.5, axis=0) l_upper2 = np.percentile(yp, 97.5, axis=0) # 做图 plt.errorbar(x,y,yerr=e, fmt="o",capsize=4) plt.plot(xi, l_fit, 'r-') plt.fill_between(xi, l_lower, l_upper, color="m",alpha=0.2) plt.fill_between(xi, l_lower2, l_upper2, color="gray",alpha=0.2) plt.plot(xi, nom - std, c='c') plt.plot(xi, nom + std, c='c') plt.plot(xi, nom - 2*std, c='c',ls="--") plt.plot(xi, nom + 2*std, c='c',ls="--") |
这样我们得到的拟合结果为 y = 1.10x -0.27 ，与不带误差拟合的结果一致，其误差为 [0.12,0.40]，与带误差拟合的结果相当。它的1σ和2σ的置信带都和 uncertainties 模块给出的边界基本吻合。
这里我们只是从数值计算的角度讨论了直线拟合中置信带的画法操作。如果对相关的数学背景和证明感兴趣，可以参考《统计推断》一书。