解决数据拟合问题最重要方法是最小二乘法和回归分析。如,我们需要从一组测定的数据(例如N个点(xi,yi)(i=0,1,…,m))去求得自变量 x 和因变量 y 的一个近似解表达式 y=f(x),这就是由给定的 N 个点(xi,yi)(i=0,1,…,m)求数据拟合的问题。(注意数据拟合和数据插值是不同的,举个例子:因为测量数据往往不可避免地带有测试误差,而插值多项式又通过所有的点(xi,yi),这样就使插值多项式保留了这些误差,从而影响逼近精度,使得插值效果不理想) 所以使用最小二乘法曲线拟合法:即寻求已知函数的一个逼近函数y=f(x),使得逼近函数从总体与已知函数的偏差按某种方法度量能达到最小,而又不一定通过全部的点(xi,yi),这个时候就需要使用最小二乘法曲线拟合法。 数据拟合的具体做法是:对给定的数据(xi,yi)(i=0,1,…,m),在取定的函数类ϕ\phiϕ中使误差r i =p(x i) −y i (i=0,1,…,m)r_{i}=p\left(x_{i}\right)-y_{i}(i=0,1, \ldots, m)ri=p(xi)−yi(i=0,1,…,m)的平方和最小,即 [∑i = 0 mri2=∑i = 0 m[ p(xi −yi )] 2] min \left[\sum_{i=0}^{m} r_{i}^{2}=\sum_{i=0}^{m}\left[p\left(x_{i}-y_{i}\right)\right]^{2}\right]_{\min }[i=0∑mri2=i=0∑m[p(xi−yi)]2]min 从几何意义讲,即寻求与给定点x i −y i (i=0,1,…,m)x_{i}-y_{i}(i=0,1, \ldots, m)xi−yi(i=0,1,…,m) 的距离平方和为最小的曲线y=p(x)。函数p(x)称为拟合函数或最小二乘解,求拟合函数p(x)的方法称为曲线拟合的最小二乘法。 在曲线拟合中,函数类ϕ\phiϕ 可有不同的选取方法。 MATLAB工具箱中提供了最小二乘拟合函数 polyfit() -->多项式曲线拟合 具体调用格式有三种:
P = polyfit(X,Y,N)[P,S] = polyfit(X,Y,N)[P,S,MU] = polyfit(X,Y,N)(1)P = polyfit(X,Y,N) 返回次数为 n 的多项式 p(x) 的系数,该阶数是 y 中数据的最佳拟合(在最小二乘方式中)。p 中的系数按降幂排列,p 的长度为 n+1. 其中 X 为输入的向量x,Y 为得到的函数值,N 表示拟合的最高次数,返回的P值为拟合的多项式:
P ( 1 ) ×X n+ P ( 2 ) ×Xn−1+ … + P ( N ) × X + P ( N + 1 ) P(1) \times X^{n}+P(2) \times X^{n-1}+\ldots+P(N) \times X+P(N+1)P(1)×Xn+P(2)×Xn−1+…+P(N)×X+P(N+1)
(2) [P,S] = polyfit(X,Y,N)还返回一个结构体 S,S可用作 polyval 的输入来获取误差估计值(S为由范德蒙矩阵的QR分解的R分量)。 其中 X 为输入的向量x,Y 为得到的函数值,N 表示拟合的最高次数,返回的P值为拟合的多项式:
P ( 1 ) ×X n+ P ( 2 ) ×Xn−1+ … + P ( N ) × X + P ( N + 1 ) P(1) \times X^{n}+P(2) \times X^{n-1}+\ldots+P(N) \times X+P(N+1)P(1)×Xn+P(2)×Xn−1+…+P(N)×X+P(N+1)
(3)[P,S,MU] = polyfit(X,Y,N)还返回 mu,mu是一个二元素向量,包含中心化值和缩放值。mu(1) 是 mean(x),mu(2) 是 std(x)。使用这些值时,polyfit 将 x 的中心置于零值处并缩放为具有单位标准差:
x ^= x−xˉσx\hat{x}=\frac{x-\bar{x}}{\sigma_{x}}x^=σxx−xˉ
这种中心化和缩放变换可同时改善多项式和拟合算法的数值属性。
下面来看一些具体的例子:(来源于帮助文档,改编)
将多项式与三角函数拟合: 在区间 [0,4*pi] 中沿余弦曲线生成 10 个等间距的点。 >> x = linspace(0,4*pi,10);>> y = cos(x);使用 polyfit 将一个 7 次多项式与这些点拟合。
>> p = polyfit(x,y,7);在更精细的网格上计算多项式并绘制结果图。
>> x1 = linspace(0,4*pi);>> y1 = polyval(p,x1);>> figure>> plot(x,y,'o')>> hold on>> plot(x1,y1,'m')>> hold off运行结果如下:
将 4 次多项式与 5 个点拟合。通常,对于 n 个点,可以拟合 n-1 次多项式以便完全通过这些点。
>> p = polyfit(x,y,4);在由 0 和 2 之间的点组成的更精细网格上计算原始函数和多项式拟合。
>> x1 = linspace(0,2);>> y1 = 1./(2+x1);>> f1 = polyval(p,x1);在更大的区间 [0,2] 中绘制函数值和多项式拟合,其中包含用于获取以圆形突出显示的多项式拟合的点。多项式拟合在原始 [0,1] 区间中的效果较好,但在该区间外部很快与拟合函数出现差异。
>> figure>> plot(x,y,'o')>> hold on>> plot(x1,y1)>> plot(x1,f1,'b--')>> legend('y','y1','f1')运行结果如下:
确定6次逼近多项式的系数。
p=polyfit(x,y,6)p =
0.0019-0.02590.1242-0.1785-0.34421.2684-0.0095为了查看拟合情况如何,在各数据点处计算多项式,并生成说明数据、拟合和误差的一个表。
f=polyval(p,x);T=table(x,y,f,y-f,'VariableNames',{'X','Y','Fit','FitError'})T =
48×4 table
XY Fit FitError _______________________________ 0 0-0.0094811 0.00948110.10.112460.11375 -0.00128760.2 0.22270.22919 -0.00649110.30.328630.33619 -0.00755990.40.428390.43431 -0.0059140.5 0.52050.52334 -0.00284080.60.603860.60326 0.000593920.7 0.67780.6742 0.00359810.8 0.74210.73643 0.00566950.90.796910.79033 0.0065794 1 0.84270.83637 0.00633161.10.880210.8751 0.00510581.20.910310.907120.0031941.30.934010.93307 0.000938261.40.952290.95361 -0.0013231.50.966110.9694 -0.0032971...在该区间中,插值与实际值非常符合。创建一个绘图,以显示在该区间以外,外插值与实际数据值如何快速偏离。
x1=(0:0.10:2*pi)'; y1=erf(x1); f1=polyval(p,x1); plot(x,y,'o') plot(x1,y1,'-') hold on plot(x1,y1,'-') plot(x1,f1,'g-') axis([0 5.5 0 3])T =
11×2 table
yearpop_____________1750 7.91e+081775 8.56e+081800 9.78e+081825 1.05e+0918501.262e+0918751.544e+091900 1.65e+0919252.532e+0919506.122e+091975 8.17e+0920001.156e+10plot(year,pop,'square') 使用带三个输入的 polyfit 拟合一个使用中心化和缩放的 5 次多项式,这将改善问题的数值属性。polyfit 将 year 中的数据以 0 为进行中心化,并缩放为具有标准差 1,这可避免在拟合计算中出现病态的 Vandermonde (范德蒙)矩阵。
使用带四个输入的 polyval,根据缩放后的年份 (year-mu(1))/mu(2) 计算 p。绘制结果对原始年份的图。
f = polyval(p,year,[],mu);hold onplot(year,f)计算在 x 中的点处拟合的多项式 p。用这些数据绘制得到的线性回归模型。
f = polyval(p,x);plot(x,y,'o',x,f,'-') legend('data','linear fit')计算以 p 为系数的一次多项式在 x 中各点处的拟合值。将误差估计结构体指定为第三个输入,以便 polyval 计算标准误差的估计值。标准误差估计值在 delta 中返回。
[y_fit,delta] = polyval(p,x,S);绘制原始数据、线性拟合和 95% 预测区间y±2Δy \pm 2 \Deltay±2Δ。
plot(x,y,'go')hold onplot(x,y_fit,'r-')plot(x,y_fit+2*delta,'m--',x,y_fit-2*delta,'m--')title('Linear Fit of Data with 95% Prediction Interval')legend('Data','Linear Fit','95% Prediction Interval')最小二乘拟合多项式系数,以向量的形式返回。p 的长度为 n+1,包含按降幂排列的多项式系数,最高幂为 n。如果 x 或 y 包含 NaN 值且 n < length(x),则 p 的所有元素均为 NaN。
使用 polyval 计算 p 在查询点处的解。
S – 误差估计结构体 (结构体) 误差估计结构体。此可选输出结构体主要用作 polyval 函数的输入,以获取误差估计值。S 包含以下字段: 字段说明RVandermonde 矩阵 x 的 QR 分解的三角因子df自由度normr残差的范数如果 y 中的数据是随机的,则 p 的估计协方差矩阵是 (Rinv*Rinv’)*normr^2/df,其中 Rinv 是 R 的逆矩阵。
如果 y 中数据的误差呈独立正态分布,并具有常量方差,则 [y,delta] = polyval(…) 可生成至少包含 50% 的预测值的误差边界。即 y ± delta 至少包含 50% 对 x 处的未来观测值的预测值。
mu – 中心化值和缩放值中心化值和缩放值,以二元素向量形式返回。mu(1) 为 mean(x),mu(2) 为 std(x)。这些值以单位标准差将 x 中的查询点的中心置于零值处。
使用 mu 作为 polyval 的第四个输入以计算 p 在缩放点 (x - mu(1))/mu(2) 处的解。
局限性
在涉及很多点的问题中,使用 polyfit 增加多项式拟合的次数并不总能得到较好的拟合。高次多项式可以在数据点之间振动,导致与数据之间的拟合较差。在这些情况下,可使用低次多项式拟合(点之间倾向于更平滑)或不同的方法,具体取决于该问题。多项式在本质上是无边界的振荡函数。所以,它们并不非常适合外插有界的数据或单调(递增或递减)的数据。算法: polyfit 使用 x 构造具有 n+1 列和 m = length(x) 行的 Vandermonde 矩阵 V 并生成线性方程组(x1nx1n−1 … 1n2x2−1 … 1 ⋮ ⋮ ⋱ ⋮nmxm−1 … 1 ) (p1p2 ⋮pn+1 )=(y1y2 ⋮ym )\left(\begin{array}{cccc} x_{1}^{n} & x_{1}^{n-1} & \ldots & 1 \\ n_{2} & x_{2}-1 & \ldots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ n_{m} & x_{m}-1 & \ldots & 1 \end{array}\right)\left(\begin{array}{c} p_{1} \\ p_{2} \\ \vdots \\ p_{n+1} \end{array}\right)=\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{array}\right)⎝⎜⎜⎜⎛x1nn2⋮nmx1n−1x2−1⋮xm−1……⋱…11⋮1⎠⎟⎟⎟⎞⎝⎜⎜⎜⎛p1p2⋮pn+1⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛y1y2⋮ym⎠⎟⎟⎟⎞
其中 polyfit 使用 p = V\y 求解。由于 Vandermonde 矩阵中的列是向量 x 的幂,因此条件数 V 对于高阶拟合来说通常较大,生成一个奇异系数矩阵。在这些情况下,中心化和缩放可改善系统的数值属性以产生更可靠的拟合。