数学建模社区-数学中国

标题: 求熟悉高斯牛顿法的大神帮忙看看我的程序 [打印本页]

作者: 和子 时间: 2016-7-12 19:20
标题: 求熟悉高斯牛顿法的大神帮忙看看我的程序
clear
%本程序用于做双高斯拟合，拟合式子为
%yi=r(1).*exp(-((x-r(2))./r(3)).^2) + r(4).*exp(-((x-r(5))./r(6)).^2);
%采用的方法是高斯-牛顿法
%x,y为做双高斯拟合的点，通过下面的式子产生
x=[0:.2:10];y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.1*randn(size(x));
%假定r初始值为1~6
r=1:6;
r=r';
y_size=size(x);
x_size=size(y);
if x_size(1)==1
x=x';
end
if y_size(1)==1
y=y';
end
yi=[];
R_square=0;
B=zeros(length(r),1);
SSE=10000000;
while 1
k=1;
%控制下系数增量的步长
for j=1:7
r1=r;
r=r+k.*B;
yi=r(1).*exp(-((x-r(2))./r(3)).^2) + r(4).*exp(-((x-r(5))./r(6)).^2);
RSSE=SSE;
yy=y-yi;
SSE=sum(yy.^2);
if RSSE>=SSE
      break;
else
      k=0.5^j;
      r=r1;
end
end
SST=sum((y-mean(y)).^2);
R_square=1-SSE/SST;
%R_square为确定系数与拟合优度有关
if R_square>0.9
   break;
end
%下面的算式是对原式做泰勒展开后省略二阶以上导数得到的，具体可参看高斯牛顿法过程
D_a1=exp(-(r(2) - x).^2./r(3).^2);
D_b1=-(r(1).*exp(-(r(2) - x).^2./r(3).^2).*(2.*r(2) - 2.*x))./r(3).^2;
D_c1=(2.*r(1).*exp(-(r(2) - x).^2./r(3).^2).*(r(2) - x).^2)./r(3).^3;
D_a2=exp(-(r(5) - x).^2./r(6).^2);
D_b2=-(r(4).*exp(-(r(5) - x).^2./r(6).^2).*(2.*r(5) - 2.*x))./r(6).^2;
D_c2=(2.*r(4).*exp(-(r(5) - x).^2./r(6).^2).*(r(5) - x).^2)./r(6).^3;
D=[D_a1 D_b1 D_c1 D_a2 D_b2 D_c2];
B=D\yy;
end

得到的结果不好，运行慢，而且很快出现
Warning: Rank deficient, rank = 1, tol =  4.079239e-17.
> In shiyan_shuanggaosi at 53
哪位大神有好的思路指点下我

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)