【例 4.9.3.3-1】 演示逆累计分布函数的应用。
v=4;xi=0.9;x_xi=chi2inv(xi,v); x=0:0.1:15;yd_c=chi2pdf(x,v); %。
plot(x,yd_c,'b'),hold on
xxf=0:0.1:x_xi;yyf=chi2pdf(xxf,v); fill([xxf,x_xi],[yyf,0],'g')
text(x_xi*1.01,0.01,num2str(x_xi))
text(10,0.16,['\\fontsize{16} x~{\\chi}^2' '(4)']) text(1.5,0.08,'\\fontname{隶书}\\fontsize{22}置信水平0.9') hold off
0.20.180.160.140.120.10.080.060.040.027.77940051015 %<8> %<9>
x~?2(4)????0.9图 4.9-5
4.10 多项式拟合和非线性最小二乘
4.10.1 多项式拟合 4.10.1.1 多项式拟合原理 4.10.1.2 拟合多项式阶数的确定 4.10.1.3 多项式拟合的MATLAB实现
【例4.10.1.3-1】实施函数拟合的较完整描述示例。
%
x=0:0.1:1;y=[2.1,2.3,2.5,2.9,3.2,3.3,3.8,4.1,4.9,5.4,5.8]; dy=0.15; for n=1:6 [a,S]=polyfit(x,y,n); A{n}=a; da=dy*sqrt(diag(inv(S.R'*S.R))); DA{n}=da'; freedom(n)=S.df; [ye,delta]=polyval(a,x,S); YE{n}=ye; D{n}=delta; chi2(n)=sum((y-ye).^2)/dy/dy; end
Q=1-chi2cdf(chi2,freedom); %
subplot(1,2,1),plot(1:6,abs(chi2-freedom),'b') xlabel('阶次'),title('chi 2与自由度')
21
subplot(1,2,2),plot(1:6,Q,'r',1:6,ones(1,6)*0.5) xlabel('阶次'),title('Q 与0.5 线') chi 2ó?×?óé?è1614121080.464200.30.20.102?×′?46002?×′?460.90.80.70.60.5Q ó?0.5 ??图 4.10-1
%
clf,plot(x,y,'b+');axis([0,1,1,6]);hold on errorbar(x,YE{3},D{3},'r');hold off title('较适当的三阶拟合')
text(0.1,5.5,['chi2=' num2str(chi2(3)) '~' int2str(freedom(3))]) text(0.1,5,['freedom=' int2str(freedom(3))]) text(0.6,1.7,['Q=' num2str(Q(3)) '~0.5']) ??êêμ±μ?èy?×?ao?65.554.543.532.521.5100.20.4Q=0.66353~0.5chi2=4.9707~7freedom=70.60.81图 4.10-2
A{3},DA{3}
ans =
0.6993 1.2005 1.8869 2.1077 ans =
1.9085 2.9081 1.2142 0.1333
4.10.2 非线性最小二乘估计 4.10.2.1 伪线性最小二乘
4.10.2.2 借助fminsearch指令进行非线性最小二乘估计
【例4.10.2.2-1】取发生信号的原始模型为y(x)?3e?0.4x?12e?3.2x。x在[0, 4]中取值；y受
到噪声0.3*(rand(n,1)-0.5) 的污染。本例演示：如何以（1）
[xydata.m]
y?a1e?a3x?a2e?a4x为模型，通
过fminsearch从受污染数据中，估计出参数a?[a(1),a(2),a(3),a(4)]?[a1,a2,a3,a4]。
22
function [x,y,STDY]=xydata(k_noise) %xydata.m x=[0:0.2:4]';
yo=3*exp(-0.4*x)+12*exp(-3.2*x); rand('seed',234)
y_noise=k_noise*(rand(size(x))-0.5); y=yo+y_noise;
STDY=std(y_noise);
[twoexps.m]
function E=twoexps(a,x,y) %twoexps.m
x=x(:);y=y(:);Y=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x); E=sum((y-Y).^2);
（2）编写M脚本文件作为主文件 [exm041022_1.m] %exm041022_1.m k_noise=0.3;
[x,y,STDY]=xydata(k_noise); a0=[1 1 1 1];
options=optimset('fminsearch'); options.TolX=0.01;
options.Display='off';
a=fminsearch(@twoexps,a0,options,x,y); chi_est=twoexps(a,x,y)/STDY^2; freedom=length(x)-length(a0); Q=1-chi2cdf(chi_est,freedom); %
y_est=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x); ych='y_e_s_t=';
a1=num2str(a(1));a2=num2str(a(2));a3=num2str(a(3));a4=num2str(a(4)); char_y_est=[ych,a1,'*exp(-',a3,'*x) + ',a2,'*exp(-',a4,'*x)']; plot(x,y,'b+');hold on,plot(x,y_est,'r');hold off,axis([0,4,0,16]) text(0.4,14,'y=3*exp(-0.4*x)+12*exp(-3.2*x)')
text(0.4,12,char_y_est),text(2.5,9,['chi2=' , num2str(chi_est)]) text(2.5,7,['freedom=' , num2str(freedom)]) text(2.5,5,['Q=' , num2str(Q)])
（3）
exm041022_1
23
16141210chi2=17.70668freedom=176Q=0.40757420y=3*exp(-0.4*x)+12*exp(-3.2*x)yest=2.8967*exp(-0.38045*x) + 11.983*exp(-3.0985*x)00.511.522.533.54图 4.10-3
【例4.10.2.2-2】以
计b?b1b2，使用fminsearch估计a?a1a2。x, y原始数据同上例。 （1）
[twoexps2.m]
function E=twoexps2(a,x,y,b) %twoexps2.m
x=x(:);y=y(:);Y=b(1)*exp(-a(1)*x)+b(2)*exp(-a(2)*x); E=sum((y-Y).^2);
（2）
[exm041022_2.m]
%exm041022_2.m k_noise=0.3;
[x,y,STDY]=xydata(k_noise); a0=[1 2]';
options=optimset('fminsearch'); options.TolX=0.001; options.Display='off'; %
while 1
Mb=exp(-x*a0'); b=Mb\\y;
a=fminsearch(@twoexps2,a0,options,x,y,b); r=norm(a-a0)/norm(a); if r<0.001;break;end a0=a; end
chi_est=twoexps2(a,x,y,b)/STDY^2; freedom=length(x)-length([a;b]); Q=1-chi2cdf(chi_est,freedom); %
y_est=b(1)*exp(-a(1)*x)+b(2)*exp(-a(2)*x); ych='y_e_s_t=';
b1=num2str(b(1));b2=num2str(b(2));a1=num2str(a(1));a2=num2str(a(2)); char_y_est=[ych ,b1,'*exp(-',a1,'*x) + ',b2,'*exp(-',a2,'*x)']; plot(x,y,'b+');hold on;plot(x,y_est,'r');hold off;axis([0,4,0,16])
??y?b1e?a1x?b2e?a2x为模型，从受污染的数据中，使用伪线性法估
?? 24
text(0.4,14,'y=3*exp(-0.4*x)+12*exp(-3.2*x)');text(0.4,12,char_y_est) text(2.5,9,['chi2=' , num2str(chi_est)]) text(2.5,7,['freedom=' , num2str(freedom)]) text(2.5,5,['Q=' , num2str(Q)])
（3）
exm041022_2 16141210chi2=18.00878freedom=176Q=0.38829420y=3*exp(-0.4*x)+12*exp(-3.2*x)yest=3.0187*exp(-0.39739*x) + 11.8737*exp(-3.1466*x)00.511.522.533.54图 4.10-4
4.10.2.3 非线性最小二乘估计指令
【例4.10.2.3-1】采用与例4.10.2.2-1相同的原始模型和受噪声污染的数据。运用lsqnonlin从受污染数据中，估计出
y?a1e?a3x?a2e?a4x的参数
a?[a(1),a(2),a(3),a(4)]?[a1,a2,a3,a4]。
（1）
[twoexps3.m]
function E=twoexps3(a,x,y)
x=x(:);y=y(:);Y=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x); E=y-Y;
（2）
clear
k_noise=0.3;
[x,y,STDY]=xydata(k_noise); a0=[1 10 0.2 1];
options=optimset('lsqnonlin');
options.TolX=0.01;options.Display='off';
a=lsqnonlin(@twoexps3,a0,[],[],options,x,y); y_est=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x);
a1=num2str(a(1));a2=num2str(a(2));a3=num2str(a(3));a4=num2str(a(4)); char_y_est=['yest=',a1,'*exp(-',a3,'*x) + ',a2,'*exp(-',a4,'*x)']; disp(['原方程 ', 'y=3*exp(-0.4*x)+12*exp(-3.2*x)']) disp(['估计方程',char_y_est])
原方程 y=3*exp(-0.4*x)+12*exp(-3.2*x)
估计方程yest=2.8381*exp(-0.37236*x) + 12.0338*exp(-3.0717*x)
4.11 插值和样条
25