initmod2.m

新的神经网络算法源程序

function [indsld, wsld] = initmod2(pmat, nbwavelon, y, nbvar);
%Initmod2: WaveNet initialization mode 2.
%
%       [indsld, wsld] = initmod2(pmat, nbwavelon, y, nbvar)
%
%Algorithm: Stepwise Regressor Selection by Orthogonalization

%By Qinghua Zhang. May, 1994.

if nbwavelon < 0  %% automatic network order determination
  nbwavelon=min(fix((size(pmat,1)-nbvar-2) / (nbvar+2)), size(pmat,2)-nbvar-1);
  AD = 1;
else AD = 0; %% network order is given
end


%initialization
resid = y'; 
STDY = sqrt(mean((y-mean(y)).^2));

% step i=1, select the bias term
li = 1;
indsld = [1];
Ij = 2:size(pmat,2);
uli = pmat(:,li)'*y';
wsld = [uli];

P = pmat(:, [1:(li-1) (li+1):size(pmat,2)]);
w_last = [pmat(:,li)];

resid = resid - uli*pmat(:,li);
NSRMSE = sqrt(mean(resid.^2)) / STDY;



if AD==1
  vary=cov(resid);
  AIC=length(y)*log(vary) + 1*4*(nbvar+2)*2;
end

for i=2:(nbwavelon+nbvar+1)
  V = pmat(:,Ij);

  P = P -  w_last * (w_last' * V);

  Pnorm2 = P.^2;
  if length(y)>1, Pnorm2=sum(Pnorm2); end;

  %If some vectors in P are nul, remove them, to avoid 0/0 later
  nulvec = find(Pnorm2 < eps);
  if ~isempty(nulvec)
    P(:,nulvec)=[]; Pnorm2(nulvec)=[]; Ij(nulvec)=[];
  end
 
  if i<=(nbvar+1) 
    litmp = 1;      % select the linear terms
  else     
    [tmpbuf, litmp] = max(((P'*y').^2)' ./ Pnorm2);
  end

  li = Ij(litmp);
  indsld = [indsld li];
  Ij = Ij([1:(litmp-1) (litmp+1):length(Ij)]);
  w_last = P(:,litmp) / sqrt(Pnorm2(litmp));

  P = P(:,[1:(litmp-1) (litmp+1):size(P,2)]);  % prepare P for next step

  uli = w_last'*y';
  wsld = [wsld uli];
  resid = resid - uli*w_last;
  NSRMSE = sqrt(mean(resid.^2)) / STDY;

  
  if i==(nbvar+1) % initialize plot   
    clf; plot([0 nbwavelon], [0 NSRMSE],'.');
    title('Stepwise Selection by Orthogonalization');
    xlabel('Nb of selected wavelets'); ylabel('NSRMSE');
    drawnow
    residhd=line('color','k','linestyle','o','erase','none', ...
        'xdata',[0],'ydata',[NSRMSE]);
    drawnow; 
  elseif i>(nbvar+1) % update plot
    set(residhd,'linestyle','*', 'xdata',[i-(nbvar+1)], 'ydata',[NSRMSE]);
    drawnow
  end

  if AD==1
    vary=cov(resid);
    N=length(y);
    dim=(nbvar+2)*(i-nbvar-1)+nbvar+1;
%dim=(0+2)*(i-nbvar-1)+nbvar+1;
    tamb=((1+dim/N)/(1-dim/N))*0.5*vary;
   AIC=[AIC tamb];
%AIC=[AIC length(y)*log(vary)+i*4*(1+2)*1.5];
  end
end

global Pause_Time; Pause_Time=0;

if AD==1 
 time_on = clock;
 AIC = AIC(nbvar+2:length(AIC));
 clf; plot(AIC);
 title('Akaike''s Final Prediction Error Criterion')
 xlabel('number of wavelets'); ylabel('FPEC');

 [tmp, nbwavelon] = min(AIC);
 nbw=input(['  Enter the number of wavelets [Default=' ...
             int2str(nbwavelon) ']: ']);

 Pause_Time = etime(clock, time_on); 

 if ~isempty(nbw)
   if(nbw<1 | nbw>length(AIC))
     disp([7 '  Invalid value. Default value is used.']);
   else   nbwavelon = nbw;
   end
 end

 indsld = indsld(1:(nbwavelon+nbvar+1));
 wsld = wsld(1:(nbwavelon+nbvar+1));
end