sre_nonline_regress.m

非线形回归程序进行最优化.(我上传的五个程序都是我个人做研究中写的,由于时间比较匆忙且只是自己写文章时做的实验,没有什么注释,如果有什么不懂的,可以email联系我:)

% nonlinear regressive anllysis

clear all
close all
%______________________________________________
%  model I: y=a* (1- exp(0-b*x) )
%         a=108; b=0.5; 
% 

% Generating exprimental data of model I.
a= 108;
b= 0.5;

% for i = 1:1:15;
%     y(i)= a* (1- exp(0- b*i) );
%     x(i)= i;
% end
% y= y + 40*abs(rand(1,15));

% load new dataset

load newkekaoxing
% 重复测试的次数
% load timebetweenfailures;
s1{1,1} = SYS1(:,1)';
s1{1,2} = SYS1(:,2)';
s1{2,1} = SS3(:,1)';
s1{2,2} = SS3(:,2)';
s1{3,1} = CSR1(:,1)';
s1{3,2} = CSR1(:,2)';
s1{4,1} = CSR2(:,1)';
s1{4,2} = CSR2(:,2)';
s1{5,1} = CSR3(:,1)';
s1{5,2} = CSR3(:,2)';
s1{6,1} = newdataset(:,1)';
s1{6,2} = newdataset(:,2)';
for KK = 1:1  %length(s1)     %  -------------------------------------------------------------------------------1 th end                            
y = s1{KK,1}; x = s1{KK,2};
end
for KK = 2:1:length(x)
    x(KK) = x(KK-1) + x(KK);
end


% estimating parameters of models
% if ppp==1   select Method A; else  select Method B.
ppp = 2;

[nry,ncy]= size(y);
  if nry<ncy   y= y';
  else   ncy= nry;  end
  
% method A: Gaussan-newton method,(reference of Fang Kaitai. p.168--171)
%*********************************************
% Step 1. initial values of estimated parameters a,and b.
fa1= [100,1];
   
fa0= [1e4,1e5];
fa0= fa0';
fa1= fa1';

ebslong= 1e-10;

f=zeros(ncy,1);  % f(i) is the estmated value of y(i)
yy= zeros(ncy,1);
J=zeros(ncy,2);  % J is differences of df/dfa0(1), df/dfa0(2), and df/dfa0(3) (i.e. dy/da, dy/db, dy/dc)  respectively.

t= 1;  % iteration times 
d(t)= sqrt(( fa1-fa0)'*(fa1-fa0) );  % stop criterion

if ppp==1     %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++     if 
    
% Step 2. iterative computing
    while d(t) > ebslong
       fa0= fa1;
       
       for i=1:1:ncy
            f(i)= fa0(1) * ( 1- exp(0- fa0(2) * x(i) ) );
          for j= 1:1:2
              if j== 1
                 J(i,j)= 1 - exp( 0- fa0(2)*x(i) );                 % dy/da
              else
               J(i,j)= fa0(1) * x(i) * exp( 0- fa0(2)*x(i) );        % dy/db
              end                          
         end
       end
       
       S(t)= (y-f)'*(y-f);   
       A= J'*J;
       A= inv(A);
       fa1= fa0+A*J'*(y-f);
       
       t= t+1;
       d(t)= sqrt((fa1-fa0)'*(fa1-fa0));
       if t>100    break;  end 
       
    end
    
    %***************************************
 else  %%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++    else
    ebslong1 = 0.1;
    %lamda=0.9;
  %*****************************************  
    while d(t) > ebslong 
       fa0= fa1;
       
       for i=1:1:ncy
            f(i)= fa0(1) * ( 1- exp(0- fa0(2) * x(i) ) );
          for j= 1:1:2
              if j== 1
                 J(i,j)= 1 - exp( 0- fa0(2)*x(i) );                 % dy/da
              else
               J(i,j)= fa0(1) * x(i) * exp( 0- fa0(2)*x(i) );        % dy/db
              end                          
         end
       end
       
       S(t)= (y-f)'*(y-f);   
       A= J'*J;
       A= inv(A);
       fa1= fa0+A*J'*(y-f);
       
       t= t+1;
       if t>100    break;  end
       
       d(t)= sqrt((fa1-fa0)'*(fa1-fa0));
       if d(t) <= ebslong & S(t-1)> 1/ebslong1
          fa1= fa1.*(1/t);   
          d(t)= sqrt((fa1-fa0)'*(fa1-fa0));   
       end             
       
   end    
   
    %************************************** 
    
end     %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ end  

%  oupput results **************************
fa = fa0
S = S;
xigma= S(length(S))/(ncy-2);
Covfa= xigma*A

%*******************************************
% caculating estimation value of y
 i=1:1:ncy;
 yyy(i)= fa0(1)* ( 1- exp( 0- fa0(2)*x(i) ) );

 % display y and yyy
 clf;
 axisy= [-1,1]; 
 axis=[0,ncy,axisy(1),axisy(2)];
 hold on;
 i= 1:1:ncy;
 plot(x(i),y(i),'-r',x(i),yyy(i),'ob');
 %i= 1:1:length(S);
 %plot(i,S(i),'or');