[问题求助] 求熟悉高斯牛顿法的大神帮忙看看我的程序

[复制链接]

字体大小: 正常放大

和子

9 主题	10 听众	132 积分

升级 16%

TA的每日心情

	奋斗 2016-7-18 14:35

签到天数: 46 天

[LV.5]常住居民I

自我介绍: GUSS

电梯直达

1^#

发表于 2016-7-12 19:20 |只看该作者 |正序浏览

|招呼Ta 关注Ta

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
哪位大神有好的思路指点下我

zan