initmod3.m

新的神经网络算法源程序

function [indsld, wsld] = initmod3(pmat, nbwavelon, y, nbvar);
%Initmod3: WaveNet initialization mode 3.
%
%       [indsld, wsld] = initmod3(pmat, nbwavelon, y, nbvar)
%
%Algorithm: Backward Elimination of Regressors

%By Qinghua Zhang. March, 1994.



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

nbcand = size(pmat,2);
nbwc = nbcand-(nbvar+1);

%initialization
resid = y';  
Ij = 1:nbcand;          CompIj=zeros(size(Ij));
indsld = [];
S_NSRMSE = [];
STDY = sqrt(mean((y-mean(y)).^2));

if AD==1, AIC=[]; end

for i=1:(nbwc-nbwavelon)
  V = pmat(:,Ij);
  if i==1 % compute M=inv(V'*V)
    M = inv(V'*V);
  else
    lastd=size(M,1);
    M11=M(1:litmp-1, 1:litmp-1);
    M12=M(1:litmp-1, litmp);
    M13=M(1:litmp-1, litmp+1:lastd);
    M21=M(litmp, 1:litmp-1);
    M22=M(litmp, litmp);
    M23=M(litmp, litmp+1:lastd);
    M31=M(litmp+1:lastd, 1:litmp-1);
    M32=M(litmp+1:lastd, litmp);
    M33=M(litmp+1:lastd, litmp+1:lastd);
    M = [M11 M13; M31 M33] - [M12; M32]*[M21 M23]/M22;
  end
  u = M * (V'*y');
  %[tmpbuf, litmp] = min((u).^2);
  [tmpbuf, litmp] = min((u(nbvar+2:length(Ij))).^2); %%## Do not select
  litmp = litmp+nbvar+1;                             %%## the linear terms
  li = Ij(litmp);             CompIj(li) = i;
  indsld = [indsld li];
  Ij = Ij([1:(litmp-1) (litmp+1):length(Ij)]);
  %wsld = [wsld uli];
  resid = y' - V*u;
  NSRMSE = sqrt(mean(resid.^2)) / STDY;
  S_NSRMSE = [S_NSRMSE NSRMSE];

  if i==1 % initialize plot
    clf; plot([nbwavelon nbwc], [0 NSRMSE],'.');
    title('Backward Elimination'); 
    xlabel('Nb of remaining wavelets'); ylabel('NSRMSE'); 
    drawnow
    residhd=line('color','k','linestyle','*','erase','none', ...
        'xdata',[nbwc],'ydata',[NSRMSE]);
  else % update plot
    axisvec = axis;
    if NSRMSE<=axisvec(4);
      set(residhd, 'xdata',[nbwc-i+1], 'ydata',[NSRMSE]);
    else
      clf; plot([nbwavelon nbwc], [0 NSRMSE],'.'); 
      title('Backward Elimination'); 
      xlabel('Nb of remaining wavelets'); ylabel('NSRMSE');
      residhd=line('color','k','linestyle','*','erase','none', ...
        'xdata',[nbwc:(-1):(nbwc-i+1)],'ydata',[S_NSRMSE]);
    end
  end
  drawnow

  if AD==1 & (nbwc-i)<maxaic
    vary=cov(resid); 
    N=length(y);
    dim=(nbvar+2)*(nbwc-i+1)+nbvar+1;
    tamb=((1+dim/N)/(1-dim/N))*0.5*vary;
    AIC=[tamb AIC]; 
  end
end

indsld = Ij;
wsld = u;

% draw the last point
V = pmat(:,Ij);
lastd=size(M,1);
M11=M(1:litmp-1, 1:litmp-1);
M12=M(1:litmp-1, litmp);
M13=M(1:litmp-1, litmp+1:lastd);
M21=M(litmp, 1:litmp-1);
M22=M(litmp, litmp);
M23=M(litmp, litmp+1:lastd);
M31=M(litmp+1:lastd, 1:litmp-1);
M32=M(litmp+1:lastd, litmp);
M33=M(litmp+1:lastd, litmp+1:lastd);
M = [M11 M13; M31 M33] - [M12; M32]*[M21 M23]/M22;
u = M * (V'*y');
resid = y' - V*u;
NSRMSE = sqrt(mean(resid.^2)) / STDY;
S_NSRMSE = [S_NSRMSE NSRMSE];
axisvec = axis;
if NSRMSE<=axisvec(4);
  set(residhd, 'xdata',[nbwc-i], 'ydata',[NSRMSE]);
else
  cla; plot([nbwavelon nbwc], [0 NSRMSE],'.');
  title('Backward Elimination'); 
  xlabel('Nb of remaining wavelets'); ylabel('NSRMSE');
  residhd=line('color','k','linestyle','*','erase','none', ...
  'xdata',[nbwc:(-1):(nbwc-i)],'ydata',[S_NSRMSE]);
end

global Pause_Time; Pause_Time=0;

if AD==1
  time_on = clock;

  vary=cov(resid);  
  N=length(y);
  dim=(nbvar+2)*(nbwc-i-1+1)+nbvar+1;      
  tamb=((1+dim/N)/(1-dim/N))*0.5*vary;
  AIC=[tamb AIC];  
  
  AIC=AIC(1:maxaic);
  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

  CompIj = CompIj .* (~(CompIj>nbwc-nbwavelon));
  indsld = find(~CompIj);
  wsld = []; %%%%%%%%%%%%
end