最小二乘法、最佳均方逼近、随机拟合及其MATLAB程序

1、2曲线拟合的线性最小二乘法及其MATLAB程序例2给出一组数据点3,y.）列入表2中，试用线性最小二乘法求拟合曲线， 估计其误差，作出拟合曲线.xi-2.53.6-1.7-1.1-0.800.11.5 2.7yi-192.9 68.04-85.50-36.15-26.52-9.10-8.43 -13.126.50解 (1)在MATLAB工作窗口输入程序 x=-2.5 -1.7-1.1-0.800.11.52.73.6;y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04;plot(x,y,r*),legend(实验数据(xi,yi

2、)xlabel(x), ylabel(y),title(数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图(略).(3)编写下列MATLAB程序计算f (x)在(七,y,)处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7-1.1-0.800.11.52.73.6;fi=a1.*x.A3+ a2.*x.A2+ a3.*x+ a4运行后屏幕显示关于ar,a2, %和a4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*

3、a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4,a4,1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4,19683/1000*a1+729/100*a2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*a3+a4编写构造误差平方和的MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/1

4、00*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4,a4,1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.A2; J=sum(fy.A2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/1

5、0)A2+(-4913/100 0*a1+289/100*a2-17/10*a3+a4+171/2)A2+(-1331/1000*a1+121/100*a2-11/10*a3+a4 + 723/2 0)+(-64/125*a1 + 16/25*a2-4/5*a3+a4 + 663/25)A2+(a4 + 91/10)+(1/1000*a1 + 1/100*a2 + 1/10*a3+a4 + 8 43/100)A2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)A2+(19683/1000 *a1+729/100*a2+27/10*a3+a4-13/2)A2+(5832/12

6、5*a1+324/25*a2 +18/5*a3+a4-1701/25)A2QJ为求a ,a ,a ,a使J达到最小，只需利用极值的必要条件 J = 0 1234Qak (k = 1,2,3,4)，得到关于a ,a ,a ,a的线性方程组，这可以由下面的MATLAB程 1234 序完成，即输入程序 syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)A2+(-4913/10 00*a1+289/100*a2-17/10*a3+a4.+171/2)A2+(-1331/1000*a1+1 21/100*a2-11/10*a3+a4+723/2

7、0)A2+(-64/125*a1+16/25*a2-4/5* a3+a4+663/25)A2+(a4+91/10)A2+(1/1000*a1+1/100*a2+1/10*a3+ a4+843/100)A2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)A2+(19683/ 1000*a1+729/100*a2+27/10*a3+a4-13/2)A2+(5832/125*a1+324/2 5*a2+18/5*a3+a4-1701/25)A2;Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3);Ja4=diff(J,a4);Ja11=s

8、imple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),运行后屏幕显示J分别对aj, a2 ,a3 ,a4的偏导数如下 Ja11= 56918107/10000*a1+32097579/25000*a2+1377283/25 00*a3+23667/250*a4-8442429/625 Ja21 = 32097579/25000*a1+1377283/2500*a2+23667/250*a3 +67*a4+767319/625 Ja31 = 1377283/2500*a1+23667/250*a2+67*a3+18/5*

9、a4-232 638/125 Ja41 =23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41=0，输入下列程序 A=56918107/10000, 32097579/25000, 1377283/2500, 23667/250; 32097579/25000,1377283/2500, 23667/250,67;1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18; B=8442429/625, -767319/625, 23263

10、8/125, -14859/25; C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下 C =5.0911-14.19056.4102-8.2574f=716503695845759/140737488355328*xA3 -7988544102557579/562949953421312*xA2 +1804307491277693/281474976710656*x -4648521160813215/562949953421312 故所求的拟合曲线为f (x) = 5.0911 x3 -14.1905 x2 + 6.4102 x - 8.2574 .(4)编写下

11、面的MATLAB程序估计其误差，并作出拟合曲线和数据的图 形.输入程序 xi=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6;y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04; n=length(xi);f=5.0911.*xi.八3-14.1905.*xi.八2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6;F=5.0911.*x.A3-14.1905.*x.A2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.-2; Ew=max(fy),

12、E1=sum(fy)/n, E2=sqrt(sum(fy2)/n) plot(xi,y,r*), hold on, plot(x,F,b-), hold off legend(数据点(xi,yi),拟合曲线 y=f(x), xlabel(x), ylabel(y),title(例2的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据3 , y )与拟合函数f的最大误差E，平均误差E和均方根误i iw1差E2及其数据点(七,七)和拟合曲线广f(0的图形(略).Ew =E1 =E2 =3.105 40.903 41.240 96函数逼近及其MATLAB程序 最佳均方逼近的MATL

13、AB主程序function yy1,a,WE=zjjfbj(f,X,Y,xx)m=size(f);n=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1);if n=length(Y)error X和Y的维数应该相同)endfor j=1:mfor k=1:mb(j,k)=0;for i=1:nb(j,k)=b(j,k)+feval(f(j,:),X(i)*feval(f(k,:),X(i);end endc(j)=0;for i=1:nc(j)=c(j)+feval(f(j,:),X(i)*Y(i);end enda=bc;WE=0;for i=1:nff=

14、0;for j=1:mff=ff+a(j)*feval(f(j,:),X(i);endWE=WE+(Y(i)-ff)*(Y(i)-ff);end if nargin=3 return; end yy=; for i=1:m l=； for j=1:length(xx) l=l,feval(f(i,:),xx(j); end yy=yy l; end yy=yy*a; yy1=yy; a=a;WE;例6. 1对数据X和Y,用函数y = 1, y = x, y = x2进行逼近，用所得到的逼 近函数计算在x = 6.5处的函数值，并估计误差.其中X=(1 3 4 5 678 9); Y=(-11-

15、13 -11 -7-1 7 17 29).解在MATLAB工作窗口输入程序 X= 13456789Y=-11-13-11-7-171729;f=fun0;fun1;fun2;yy,a,WE=zjjfbj(f,X,Y,6.5) 运行后屏幕显示如下 yy = 2.75000000000003 a = -7.00000000000010-4.999999999999951.00000000000000 WE = 7.172323350269439e-027 例 6.2 对数据 X 和 Y，用函数 y = 1, y = x, y = x2，y = cosx，y = ex， y = sinx进行逼近，其

16、中 X=(0 0.50 1.00 1.50 2.00 2.50 3.00)，Y二(0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645).解 在MATLAB工作窗口输入程序 X= 00.501.001.502.002.503.00;Y=00.47940.84150.98150.91260.59850.1645; f=fun0;fun1;fun2;fun3;fun4;fun5;xx=0: 0.2:3; yy,a,WE=zjjfbj(f,X,Y,xx),plot(X,Y,ro,xx,yy,b-) 运行后屏幕显示如下(图略) yy = Columns 1 throu

17、gh 7 -0.00050.20370.39390.56560.71410.83480.9236Columns 8 through 14 0.97710.99260.96910.90690.80800.67660.5191Columns 15 through 16 0.34440.1642a = 0.38280.4070-0.39010.0765-0.45980.5653WE = 1.5769e-004即，最佳逼近函数为y=0.3828+0.4070*x-0.3901*x2+0.0765*exp(x)-0.4598*cos(x)+0.5653*sin(x).8随机数据点上的二元拟合及其MATL

18、AB程序例8设节点(X,Y,Z)中的X和Y分别是在区间-3, 3和-2.5, 3.5上的50 个随机数，Z是函数Z=7-3x3e -x2-y2在(X,Y)的值，拟合点(X”Y?中的 XI=-3:0.2:3, YI=-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三 角基线性内插法、三角基三次内插法和MATLAB 4网格化坐标方法计算在(X”Y? 处的值，作出它们的图形，并与被拟和曲面进行比较.解 (1)最近邻内插法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x，y .X=-3+(3-(-3)*x;%利用x生成的随机变量.

19、Y=-2.5+(3.5-(-2.5)*y; %利用y生成的随机变量.Z = 7-3* X.A3 .* exp(-X.A2 - Y.A2); % 在每个随机点(X,Y) 处计算Z的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1); % 将坐标(XI,YI)网格化.ZI=griddata(X,Y,Z,XI,YI, nearest ) %计算在每个插值点 (XI,YI)处的插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z),title(用最近邻内插法拟合函数z =7-3 x

20、A3 exp(-xA2 - yA2) 的曲面和节点的图形)%legend(拟合曲面，节点 (xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示用最近邻内插法拟合函数Z=7-3x3e-x2-y2在两组不同节点处的曲 面及其插值4 (略).三角基线性内插法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x，y .X=-3+(3-(-3)*x;%利用x生成 上的随机变量.Y=-2.5+(3.5-(-2.5

21、)*y; %利用y生成 上的随机变量.Z = 7-3* X.A3 .* exp(-X.A2 - Y.A2); % 在每个随机点(X,Y) 处计算Z的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1);% 将坐标(XI,YI )网格化.ZI=griddata(X,Y,Z,XI,YI, linear ) %计算在每个插值点 (XI,YI )处的插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z), title(用三角基线性内插法拟合函数z =7-3 xA3 exp(-xA2

22、- yA2)的曲面和节点的图形)%legend(拟合曲面，节点 (xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示用三角基线性内插法拟合函数Z=7-3x3e -x2y在两组不同节点处 的曲面和节点的图形及其插值4 (略).三角基三次内插法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x，y .X=-3+(3-(-3)*x;%利用x生成 上的随机变量.Y=-2.5+(3.5-(-2.5)*y; %

23、利用y生成 上的随机变量.Z = 7-3* X.A3 .* exp(-X.A2 - Y.A2); % 在每个随机点(X,Y) 处计算Z的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1);% 将坐标(XI,YI)网格化.ZI=griddata(X,Y,Z,XI,YI, cubic ) %计算在每个插值点 (XI,YI)处的插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z), title(用三角基三次内插法拟合函数z =7-3 xA3 exp(-xA2 - yA2)的曲面

24、和节点的图形)%legend(拟合曲面，节点 (xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示用三角基三次内插法拟合函数Z=7-3x3e -x2-y2在两组不同节点处 的曲面和节点的图形及其插值4 (略).MATLAB 4网格化坐标方法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x， y .X=-3+(3-(-3)*x;%利用x生成 上的随机变量.Y=-2.5+(3.5-(-2.5)*y; %利用y生成 上的随机变量.Z = 7-3* X.A3 .* exp(-X.A2 - Y.A2); % 在每个随机点(X,Y) 处计算Z的值.X1=-3:0.2:3; Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1);% 将坐标(XI,YI )网格化.ZI=griddata(X,Y,Z,XI,YI, v4 ) %计算在每个插值点 (XI,YI)处的插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), yl