基于最小二乘法的多项式拟合

基于最小二乘法的多项式拟合 一 最小二乘法的基本原理 y,f(x)设已知某物理过程的一组观测数据 (x,f(x))i,1,2,？,mii ， . (1) x,(x)y,f(x),(x)i要求在某特定函数类寻求一个函数作为的近似函数，使得二者在上的残差 ,,,(x),f(x)i,1,2,？,miii， (2) 按某种度量标准为最小，这就是拟合问题. T,,[,,,,？,,],,01iim要求残差按某种度量标准为最小，即要求由残差构成的残差向量的某种范 ,,,1,数为最小，要求，或即 mm ,,,,,(x),f(x),,iii1,,00ii ,,max,,max,(x),f(x)iii,ii 为最小，这本来都是很自然的，可是计算不太方便.通常要求: 11mm2222,,(,),{[,(x),f(x)]},,iii2,0,0ii 或者 mm222,,,,[,(x),f(x)],,iii2,,00ii (3) 为最小.这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法.就是说，最小二乘法提供了一种数学方 ,法，利用这种可以对实验数据实现在最小平方误差意义下的最好拟合.在曲线拟合中，函数类可有不同的选取方法.下面就常用的多项式拟合做介绍。 二 多项式拟合 (x,y)n(n,m)ii,假设给定数据点(i=0,1,…,m)，为所有次数不超过的多项式构成的函数类，现求一 nkp(x),ax,,,nkk,0,使得 2mmn,,2k,,I,p(x),y,ax,y,min,,,,,niikiii,,,00ik0,, (4) p(x)n当拟合函数为多项式时，称为多项式拟合，满足式(4)的称为最小二乘拟合多项式。特别地，当n=1时，称为线性拟合或直线拟合。 显然 mn2kI,(ax,y),,kii,,00ik a,a,？aI,I(a,a,？a)01n01n为的多元函数，因此上述问题即为求的极值 问题。由多元函数求极值的必要条件，得 mn,Ikj,2(ax,y)x,0,j,0,1,？,n,,kiii,a,,00ikj (5) 即 nmm，jkj(x)a,xy,j,0,1,？,n,,,ikii,,00,0kii (6) a,a,？a01n(6)是关于的线性方程组，用矩阵表示为 mmm，，，，n，mxx1？y,,ii,,,i,,ai,0i,0，，i,00,,,,mmmm,,2n，1,,,,axxx？1xy,,,iii,,,ii,,,,,i,0i,0i,0i,0,,？,,,,？？？？,,,,m,,mmman，，nnn，12n,,,,xyxxx？,ii,,,iii,,,,i,0，，i,0i,0i,0，， (7) 式(6)或式(7)称为正规方程组或法方程组。 ak可以证明，方程组(7)的系数矩阵是一个对称正定矩阵，故存在唯一解。从式(7)中解出(k=0,1,…，n)，从而可得多项式 nkp(x),ax,nk,0k (8) m2,,p(x),y,niip(x)p(x)nn,0i可以证明，式(8)中的满足式(4)，即为所求的拟合多项式。我们把称 p(x)n为最小二乘拟合多项式的平方误差，记作 m22,,r,p(x),y,nii2,0i 由式(5)可得 mnm22kr,y,a(xy),,,ikii2,,,000iki (9) 多项式拟合的一般方法可归纳为以下几步: (1) 由已知数据画出函数粗略的图形——散点图，确定拟合多项式的次数n; mmjjx(j,0,1,？,2n)xy(j,0,1,？,2n),,iii,0,0ii(2) 列表计算和; a,a,？a01n(3) 写出正规方程组，求出; nkp(x),ax,nk,0k(4) 写出拟合多项式。 n,mn,mn,m在实际应用中，或;当时所得的拟合多项式就是拉格朗日或牛顿插值多项式。 TR(,)ii例1 测得铜导线在温度(?)时的电阻如表6-1，求电阻R与温度 T的近似函数关系。 i 0 1 2 3 4 5 6 19.1 25.0 30.1 36.0 40.0 45.1 50.0 Ti(?) 76.30 77.80 79.25 80.80 82.35 83.90 85.10 R(,)i 解 画出散点图(图6-2)，可见测得的数据接近一条直线，故取n=1，拟合函数为 R,a，aT01 列表如下 i 2TRTRTiiiii 0 19.1 76.30 364.81 1457.330 1 25.0 77.80 625.00 1945.000 2 30.1 79.25 906.01 2385.425 3 36.0 80.80 1296.00 2908.800 4 40.0 82.35 1600.00 3294.000 5 45.1 83.90 2034.01 3783.890 6 50.0 85.10 2500.00 4255.000 245.3 565.5 9325.83 20029.445 , 正规方程组为 a7245.3565.5，，，，，，0,,,,,,,a245.39325.8320029.4451，，，，，， 解方程组得 a,70.572,a,0.92101 故得R与T的拟合直线为 R,70.572，0.921T 利用上述关系式，可以预测不同温度时铜导线的电阻值。例如，由R=0得T=-242.5，即预测温度 T=-242.5?时，铜导线无电阻。 6-2 例2 已知实验数据如下表 i 0 1 2 3 4 5 6 7 8 1 3 4 5 6 7 8 9 10 xi 10 5 4 2 1 1 2 3 4 yi 试用最小二乘法求它的二次拟合多项式。 解 设拟合曲线方程为 2y,a，ax，ax012 列表如下 I 2342xyxyxxxxyiiiiiiiii 0 1 10 1 1 1 10 10 1 3 5 9 27 81 15 45 2 4 4 16 64 256 16 64 3 5 2 25 125 625 10 50 4 6 1 36 216 1296 6 36 5 7 1 49 343 2401 7 49 6 8 2 64 512 4096 16 128 7 9 3 81 729 6561 27 243 8 10 4 100 1000 10000 40 400 53 32 381 3017 25317 147 1025 , 得正规方程组 952381a32，，，，，，0,,,,,,523813017a,1471,,,,,, ,,,,,,381301725317a10252，，，，，， 解得 a,13.4597,a,,3.6053a,0.2676012 故拟合多项式为 2y,13.4597,3.6053，0.2676x 三 多项式拟合的MATLAB实现 用Polyfit函数P=polyfit(x,y,n)对数据进行拟合，返回n次多项式的系数，并用降序排列的向量表示，长度为n+1. nn,1 p(x),px，px，？，px，p12nn，1 [p,s]=polyfit(x,y,n)返回多项式系数向量p和矩阵s。s与polyval函数一起用时，可以得到预测值的误差估计。 例如在MATLAB界面中输入一下命令 >> x=[0 0.0385 0.0963 0.1925 0.2888 0.385]; >> y=[0.042 0.104 0.186 0.338 0.479 0.612]; >> [p,s,mu]=polyfit(x,y,5) 输出结果为: p = Columns 1 through 5 0.0193 -0.0110 -0.0430 0.0073 0.2449 Column 6 0.2961 54320.0193x,0.0110x,0.043x，0.0073x，0.2449x，0.2961说明拟合的多项式为: s = R: [6x6 double] df: 0 normr: 2.3684e-016 mu = 0.1669 0.1499 自由度为 0 标准偏差为 2.3684e-016