最小二乘法的基本原理和多项式拟合matlab实现

1、最小二乘法的基本原理和多项式拟合从整体上考虑近似函数P(x)同所给数据点(Xi，yD(i=0,1,m)误差i=P(Xi)-yi(i=0,i, . ,m)ri= p(xi)yi(i=0,1,m)绝对值的最大值maxmri，即误差 向量 mr =(r0,i,rm)T的8范数；二是误差绝对值的和 匕ri，即误差向量r的1 一 my r 2范数；三是误差平方和 菖i的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便丁微分运算，后一种方法相当丁考虑2范数的平方, m因此在曲线拟合中常采用误差平方和 乙来 度量误差ri(i=0 , 1,，m)的整数据拟合的具体作法是：对给定数据(xi,yi)

2、 (i=0,1,，m),在取定的函数类中中,求p(x短小，使误差ri=p(xi)yi(i=0,1, . ,m)的平方和最小，即'U2p(xi) - yi 2 =mini -0i =0从几何意义上讲，就是寻求与给定点(xi,yi)(i=0,1,m)的距离平方和为最y =p(x)(图6-1)。函数p(x)称为拟合函数或最小二乘解，求拟合 函数p(x)的方法称为曲线拟合的最小二乘法。中可有不同的选取方法.6 1假设给定数据点（xi，%）（i=0,1,m），中为所有次数不超过n（nm）的多项式构nkPn（x） = akX成的函数类，现求一心，使得mI - ' Ipn(Xi) - yii

3、 £=£ fz akX： yj =min i =0 >=0当拟合函数为多项式时，称为多项式拟合，满足式（1）的Pn（X）称为最小二乘拟合多项式。特别地，当n=1m nI = L (' a：x： - yD2i z0 k z0为瓦,礼色的多元函数，因此上述问题即为求1 =l（a°,a1,an ）的极值问题。由多元函数求极值的必要条件，得m n=2、(' akXik -yi)xj =0, a j ik =0j =0,1,nn m('、' Xi k=0i =0m j "kji)ak =Xi yi ,i =0j =0,1, ,

4、n(3)是关丁 a0,a1,色的线性方程组，用矩阵表示为m、Xii Wm' Xii =0 m'、Xi2i Wm一 n I£ Xii=0mn书z Xii =0；1£ Vi=0m£ xyii =0m n 1Xii Wmx - 2nXii =0m寸 nx yij=0式（3）或式（4可以证明，方程组（4）的系数矩阵是一个对称正定矩阵,故存在唯一解。从式（4）中解出ak（k=0,1,，n）nPn(x) = .，akxk k=0可以证明，式(5)中的Pn(X)满足式(1),即Pn(x)为所求的拟合多项式。我mw bn(xi) - yi f 们把i A称为最小二

5、乘拟合多项式Pn(X)r|2 =Z hn(xi) yi z0|2由式2 mnmr|2 = L y：ak(' x")i =0k z0i z0多项式拟合的一般方法可归纳为以下几步:由已知数据画出函数粗略的图形散点图,确定拟合多项式的次数n；(j =0,1,2n)Z xij (j =0,1，,2n) z xijyi歹0表计算仕和i=0 写出正规方程组，求出ao，a1"annPn(x) =akXk写出拟合多项式k&在实际应用中，n<m或n<m ；当n = m时所得的拟合多项式就是拉格朗日或牛顿插值多项式。例1测得铜导线在温度Ti(C )时的电阻Ri(Q)

6、如表6-1 ,求电阻R与温度Ti0123456T(C)19.125.030.136.040.045.150.0Ri (。)76.3077.8079.2580.8082.3583.9085.10解 画出散点图(图6-2)，可见测得的数据接近一条直线，故取 n=1,拟合函 数为R = a° ' aT列表如下iTiRiTi2TiR019.176.30364.811457.330125.077.80625.001945.000230.179.25906.012385.425336.080.801296.002908.800440.082.351600.003294.000545.18

7、3.902034.013783.890650.085.102500.004255.000E245.3565.59325.8320029.4457245.3a0565.5245.3 9325.83 涌1120029.445ao = 70.572,a = 0.921故得R与TR =70.572 0.921T利用上述关系式，可以预测不同温度时铜导线的电阻值。例如，由 R=0得 T=-242.5T=-242.5 80 -.J1 0030T6-2例2例2i012345678Xi1345678910y1054211234解设拟合曲线方程为2y =a° +ax + a2XIXiV2Xi3Xi4Xi

8、xyi2Xi yi01101111010135927811545244166425616643522512562510504613621612966365714934324017496826451240961612879381729656127243810410010001000040400E53323813017253171471025a。=13"952381a0523813017 a1381 3017 25317La2 一>.4597, a = -3.6053一 32 1147J025 一a2 -0.26762y =13.4597-3.6053 0.2676x2*三定理1

9、设节点X0，Xl"，Xn互异，则法方程组（4证 由克莱姆法则，只需证明方程组（4用反证法，设方程组（4m'、'Xii -0m'、'Xi2i有非零解。将式(8)m- nXii J0n mk =0i =0m' Xin1i £X，)ak =0,-a。aaZ Vi=0 mZ Xi yii=0aan_mn1I Xi 乂J=0_mm- 2n' Xiimn' Xii或mn 1'Xii £j =0,1, ,n(7(8中第j个方程乘以aj(j=0,1,，n),然后将新得到的n+1个方程左n n m'、a广 0右

10、两端分别相加，得H YXij k)ak0 =0k=0 i=0nn mm n nm nnmZ ajz (E Xij*)ak =Z Z Z akajX：* =£ (Z ajXij)(E akXk)=ZPn(X)j =0 |Lk=0i=0i=0j=0kz0i =0 j =0k=0i=0nPn(X) _akXkk -0pn(Xi) = 0 (i=0,1,m)pn(X)是次数不超过n的多项式，它有m+bn个相异零点，由代数基本定理，必须有a0 = a1an = 0 ,与齐次方程组有非零解的假设矛盾。因此正规方程组(4)n.Pn(X) = ' akX必有唯一解。定理2 设a0,a1,an

11、是正规方程组(4)的解，贝UkT是满足式(1n证只需证明，对任意一组数b0,b1,，bn组成的多项式Qn(X)=匕bkX ,何有mm"Qn(Xi) -yiPn(Xi)-yii =0i =0mm'、Qn (xi ) yi 1 - ' ，Pn (xi ) y i zQi zQmm=、Qn(Xi)"n(Xi)2 2、 Qn(Xi) - Pn(Xi)】-Pn(Xi) - 乂 1i zQi zQm n_ n"I nm,n云。+2££ (bj aDxJakX yi |=密(加aj £Z akX _ yi 时 *ig j=QJ=Qj

12、=Q、i=0 -kJ Jj因为ak(k=0,1,，n)是正规方程组(4)的解，所以满足式(2mm'、Qn(Xi) -yjIpn(Xi) - yi_0i =0i=Q故 Pn(X)*四在多项式拟合中，当拟合多项式的次数较高时，其正规方程组往往是病态的。而 拟合节点分布的区间'X0，Xm Xi(i=0,1,，m)不使用原始节点作拟合，将节点分布区问作平移，使新的节点Xi关丁原点对Xi = Xi - X。Xm ,i =0,1, ,m2 对平移后的节点Xi(i=0,1,，m),X = pXi , i =0,1, ,mp =2r(m 1) 3)2r(r是拟合次数)(10(11经过这样调整可

13、以使Xi的数量级不太大也不太小，特别对丁等距节点X =X。+ih （i =0,1，m），作式（10）和式（11）两项变换后，其正规方程 组的系数矩阵设 为A,则对14次多项式拟合，条件数都不太大，都可以得到拟合次数1234cond2 (A)=1<9.9<50.3<435在实际应用中还可以利用正交多项式求拟合多项式。一种方法是构造离散正交 多项式；另一种方法是利用切比雪夫节点求出函数值后再使用正交多项式。这两种方法都使正规方程 组的系数矩阵为对角矩阵，从而避免了正规方程组的病态。例如 m=19, xo =328,h=1, x1 = xo+ih , i=0,1,，19,即节点 分

14、布在328,347, 直接用Xi构造正规方程组系数矩阵Acond2(A。)=2.25 1016328 347i =0,1, ,19为=xi -2X构造正规方程组系数矩阵Acond2(A) =4.483868 1016比 concHA)降低了 13P =4200.149819")4i =0Xi = pxi ,i =0,1, ,19用X 构造正规方程组系数矩阵A2cond2(A2) =6.839乂比 cond2(A1)降低了 3如有必要，在得到的拟合多项式 Pn(xD中使用原来节点所对应的变量x,可写为Qn(X)= Pn(P (X _Xo 2Xm)仍为一个关丁 x的n次多项式，正是我们要

15、求的拟合多项式。Matlab 实现：在MATLAB中可用函数aa=ployfit(x,y,n) 来求得参数a，从而实现线性 最小二乘拟合(当然也可以实现多项式函数拟合)，其中参数x为给定结点数据 的横坐标向量，y为对应的纵坐标向量，n为多项式的次数(如线性拟合则n为 1,二次多项式拟合为2等)，函数返回值aa为拟合的多项式系数。在求得多项式系数后，为了求得多项式的值，可用MATLAB函数y=polyval(aa,x) 求得系数为aa的多项式在指定点x的函数值y。(1)【例】用多项式最小二乘拟合求电阻R与温度t之间的关系R = at + b :%多项式最小二乘数拟合(线性拟合实例)>>

16、;t=20.5 32.5 51 73 95.7?r=765 826 873 942 1032?>>aa=polyfit(t,r,1) ?%插值点，其中1表示一次多项式，即直线>>a=aa(1) %插值求得一次项系数a=3.3940>>b=aa(2) %插值求得常数项系数b=702.4918>>y=polyval(aa,t)%求得拟合出来的多项式在输入值为x下的y值>>plot(t,r, ' k+' , t,y, ' r')图1-1(2)例题对以下数据分别作二次和三次多项式拟合，求得多项式系数，并画 出图

17、形。X24567S10121416y648,49.289.59.79.3610,210.410.510.6二次多项式拟合程序如下：(程序中如果想显示结果就不加分号，图 1-2 ) %多项式最小二乘法拟合，参照(matlab实验实验指导书李新平' 实验六) 自己做的%多项式表示为y=ax.A2+bx+cx=2 4 5 6 7 8 10 12 14 16; r=6.4 8.4 9.28 9.5 9.7 9.86 10.2 10.4 10.5 10.6;xs=polyfit(x,r,2); %求的系数矩阵这里是向量xs,其中2表示二次多项式a=xs(1)%插值求得二次项系数b=xs(2)%插

18、值求得一次项系数c=xs(3)%插值求得常数项y=polyval(xs,x) %求得拟合出来的多项式在输入值为 x下的y值plot(x,r,'r*' ,x,y,'b-')%画出通过数据点而拟合出来的曲线三次多项式拟合程序如下：(图1-3)%多项式最小二乘法拟合，参照(matlab实验实验指导书李新平'实验六)自己做的% 多项式表示为 y=ax.A3+bx.A2+cx+dx=2 4 5 6 7 8 10 12 14 16; r=6.4 8.4 9.28 9.5 9.7 9.86 10.2 10.4 10.5 10.6;xs=polyfit(x,r,3);a=xs(1)b=xs(2)c=xs(3)d=xs y=polyval(xs,x)