最小二乘法、最佳均方逼近、随机拟合及其MATLAB程序
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
2曲线拟合的线性最小二乘法及其MATLAB?序
例2给出一组数据点(X i, y i)列入表2中,试用线性最小二乘法求拟合曲线, 估计其误差,作出拟合曲线•
解 (1)在MATLAB工作窗口输入程序
>> x=[-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];
plot(x,y , 'r*'),
lege nd( '实验数据(xi,yi)' )
xlabel( 'x' ), ylabel( 'y'),
title( '数据点(xi,yi) 的散点图’)
运行后屏幕显示数据的散点图(略)
(3)编写下列MATLAB程序计算f(x)在(X j,yj处的函数值,即输入程序
>> syms al a2 a3 a4
x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];
fi=a1.*x.A3+ a2.*x.A2+ a3.*x+ a4
运行后屏幕显示关于a1,a2, a3和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*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.50 68.04];
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*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/10)A2+(-4913/100
0*a1+289/100*a2-17/10*a3+a4+171/2)A2+(-1331/1000*a1+121/1
00*a2-11/10*a3+a4+723/20F2+(-64/125*a1+16/25*a2-4/5*a3+a
4+663/25F2+(a4+91/10F2+(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/2F2+(5832/125*a1+324/25*a2 +18/5*a3+a4-1701/25F2
为求31,32,33,34使J达到最小,只需利用极值的必要条件-丄 0 a k (k 1,2,3,4),得到关于31,32,33,34的线性方程组,这可以由下面的MATLAB程序完成,即输入程序
>> syms a1 a2 a3 a4 J=(-125/8*a1+25/4*32-5/2*a3+34+1929/10)A2+(-4913/10
00*a1+289/100*32-17/10*33+34...+171/2)A2+(-1331/1000*a1+1 21/100*a2-
11/10*a3+34+723/20)A2+(-64/125*31+16/25*32-4/5*
a3+34+663/25)A2+(34+91/10)A2+(1/1000*31+1/100*32+1/10*a3+
a4+843/100)A2+(27/8*31+9/4*32+3/2*a3+34+328/25)A2+(19683/ 1000*a1+729/100*32+27/10*a3+34-
13/2)A2+(5832/125*31+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=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3),
Ja41=simple(Ja4),
运行后屏幕显示J分别对31, 32 ,33 ,34的偏导数如下
Ja1仁
56918107/10000*31+32097579/25000*32+1377283/25
00*33+23667/250*34-8442429/625
J321 =
32097579/25000*31+1377283/2500*32+23667/250*33 +67*34+767319/625 J331 = 1377283/2500*31+23667/250*32+67*33+18/5*34-232
638/125
J341 = 23667/250*31+67*32+18/5*33+18*34+14859/25
解线性方程组J311 =0,J321 =0,J331 =0,J341 =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, 232638/125, -14859/25];
C=B/A, f=poly2sym(C)
运行后屏幕显示拟合函数f及其系数C如下
C = 5.0911 -14.1905 6.4102 -8.2574
f=716503695845759/140737488355328*xA3
-7988544102557579/562949953421312*xA2
+1804307491277693/281474976710656*x -4648521160813215/562949953421312 故所求的拟合曲线为
f (x) 5.0911 x314.1905 x2 6.4102 x 8.2574 .
(4)编写下面的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.A3-14.1905.*xi.A2+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.A2; Ew=max(fy),
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)的图形')
运行后屏幕显示数据(X i,yj与拟合函数f的最大误差平均误差E i和均方根误差E2及其数据点(X j,yj和拟合曲线y=f(x)的图形(略).