MATLAB 多项式拟合和非线性最小二乘

作者：未知信息来源:未知 2006-1-26

字体大小:小 中 大网友评论条 进入论坛

10多项式拟合和非线性最小二乘5。1 多项式拟合5。3 多项式拟合的MATLAB实现【 * 例 5。3-1 】实施函数拟合的较完整描述示例。...

5.10多项式拟合和非线性最小二乘

5.10.1 多项式拟合

5.10.1.3 多项式拟合的MATLAB实现

【 * 例 5.10.1.3-1 】实施函数拟合的较完整描述示例。
在进行函数拟合时，有几个问题要回答：
（1）采用什么函数模型？
（2）模型的结构数什么？
（3）参数的估计值如何计算？
（4）估计参数的离差？

本例采用多项式模型，模型阶数（结构参数）通过量确定，多项式系数运用最小二乘估计，并给出相应的离差。具体如下。
% 被拟合的原始数据
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; % 原数据 y 的标准差
for n=1:6 % 依次用 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 与自由度 ')
subplot(1,2,2),plot(1:6,Q,'r',1:6,ones(1,6)*0.5)
xlabel(' 阶次 '),title('Q 与 0.5 线 ')

图 5.10.1.3-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'])

图 5.10.1.3-2 三阶多项式的拟合情况

A{3},DA{3} % 拟合多项式的系数和离差
ans =
0.6993 1.2005 1.8869 2.1077
ans =
1.9085 2.9081 1.2142 0.1333

〖说明〗由于原始数据较少，被估多项式系数有较大的离差。

5.10.2 非线性最小二乘估计

5.10.2.2 借助fminsearch指令进行非线性最小二乘估计
【 * 例 5.10.2.2-1 】取发生信号的原始模型为。在 [0, 4] 中；受到噪声 0.3*(rand(n,1)-0.5) 的污染。本例演示：如何以为模型，通过 fminsearch 从受污染数据中，估计出参数。

（1）编写两个 M 函数文件：仿真数据产生文件和二乘残差计算文件。
[xydata.m]
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 脚本文件作为主文件
[exm051022_1.m]
%exm051022_1.m 利用 fminsearch 指令进行非线性参数估计
k_noise=0.3; % 控制噪声水平
[x,y,STDY]=xydata(k_noise); % 运行仿真数据产生程序，产生数据
a0=[1 1 1 1]; % 被估参数的初试猜测
options.TolX=0.01; % 控制被估参数的迭代精度
options.Display= 'off' ; % 避免显示收敛信息
a=fminsearch( 'twoexps' ,a0,options,x,y); % 计算二乘残差最小时的参数估计
chi_est=twoexps(a,x,y)/STDY^2; % 估计参数下的 Chi2 量计算
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=' ; % 注意：该格式使图形中 est 三个字母成为下标
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）以上文件都存放在用户自己的工作目录上，再使该目录成为当前目录，则在 Notebook 中就可以进行以下指令的运行，并得到如图 5.10.2.2-1 所示的图形结果。

exm051022_1

图 5.10.2.2-1 借助 fminsearch 指令估计非线性模型参数

【 * 例 5.10.2.2-2 】以为模型，从受污染的数据中，使用伪线性法估计，使用 fminsearch 估计。

（1）编写两个 M 函数文件：仿真数据产生文件（与上例相同）和二乘残差计算文件。
[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）编写 M 脚本文件作为主文件
[ exm051022_2.m ]
% exm051022_2 .m 利用伪线性和 fminsearch 混合法进行非线性参数估计
k_noise=0.3; % 控制噪声水平
[x,y,STDY]=xydata(k_noise); % 运行仿真数据产生程序，产生数据
a0=[1 2]'; % 非线性参数的初试猜测
options.TolX=0.001; % 控制被估参数的迭代精度
options.Display= 'off' ; % 避免显示收敛信息
% 下面是伪线性和 fminsearch 交互使用的迭代循环
while 1
Mb=exp(-x*a0'); % 利用初试猜测 a0 ，构成线性估计的系数矩阵。
b=Mb\y; % 线性最小二乘求线性参数 b
a=fminsearch( 'twoexps2' ,a0,options,x,y,b); % 利用 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; % 估计参数下的 Chi2 量计算
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=' ; % 注意：该格式使图形中 est 三个字母成为下标
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])
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）以上文件都存放在用户自己的工作目录上，再使该目录成为当前目录，则在 Notebook 中就可以进行以下指令的运行，并得到如图 5.10.2.2-2 所示的图形结果。