pcrcvbg.m

偏最小二乘算法在MATLAB中的实现

function [press,cumpress,minpc,b,p,qlim,t2lim,tvar] = pcrcvbg(x,y,split,pc,np,mc,gui,action)
%PCRCVBG Cross validation for PCR with contiguous blocks for use by MODLGUI

% Copyright
% Eigenvector Technologies
% 1995
global cumpress press minpc b p qlim t2lim t ssq eigs bbr tvar
handl   = get(gui,'UserData');
[mx,nx] = size(x);
[my,ny] = size(y);
if strcmp(action,'crval')
  if mx ~= my
    error('Number of samples must be the same in both blocks')
  end
  if nargin < 5
    np = 1;
  end
  if nargin < 6
    mc = 1;
  end
  press = zeros(split*ny,pc);
  ind = ones(split,2);
  for i = 1:split
    ind(i,2) = round(i*mx/split);
  end 
  for i = 1:split-1
    ind(i+1,1) = ind(i,2) +1;
  end
  for i = 1:split
    s = sprintf('Now working on test set %g out of %g',i,split);
	set(handl(2,4),'String',s)
    if mc ~= 0
      [calx,mnsx] = mncn([x(1:ind(i,1)-1,:); x(ind(i,2)+1:mx,:)]);
      testx = scale(x(ind(i,1):ind(i,2),:),mnsx);
      [caly,mnsy] = mncn([y(1:ind(i,1)-1,:); y(ind(i,2)+1:mx,:)]);
      testy = scale(y(ind(i,1):ind(i,2),:),mnsy);
    else
      calx = [x(1:ind(i,1)-1,:); x(ind(i,2)+1:mx,:)];
      testx = x(ind(i,1):ind(i,2),:);
      caly = [y(1:ind(i,1)-1,:); y(ind(i,2)+1:mx,:)];
      testy = y(ind(i,1):ind(i,2),:);
    end
    [t,p,bbr] = pcrg(calx,caly,pc);
    for j = 1:pc
      ypred = testx*bbr((j-1)*ny+1:j*ny,:)';
	  press((i-1)*ny+1:i*ny,j) = sum((ypred-testy).^2)';
    end
	figure(handl(1,1))
    plot(press((i-1)*ny+1:i*ny,:)')
    s = sprintf('PRESS for Test Set Number %g out of %g',i,split);
    title(s)
    xlabel('Number of Principal Components')
    ylabel('PRESS')
    drawnow
  end
  pause(2)
  cumpress = sum(press);
  figure(handl(1,1))
  plot(cumpress)
  title('Cumulative PRESS as a Function of Number of Principal Components')
  xlabel('Number of Principal Components')
  ylabel('PRESS')
  drawnow
  [a,minpc] = min(cumpress);
  s = sprintf('Minimum Cumulative PRESS is at %g pcs',minpc);
  figure(gui)
  s = sprintf('Minimum Cumulative PRESS is at %g pcs',minpc);
  if minpc>1
    s1 = ' Click slider to choose different number of LVs - then RESUME to continue';
    set(handl(3,6),'Max',minpc,'Value',minpc);
  else
    s1 = ' - Click RESUME to continue';
	set(handl(3,6),'Value',minpc);
	set(handl([3 5 6],6),'Visible','Off')
  end
  set(handl(4,6),'String',num2str(minpc))
  figure(gui)
  set(handl(2,4),'String',[s,s1])
  set(handl(4,12),'Visible','On')
  set(handl(2,6),'String','Latent Variables');
  set(handl(6,6),'String',num2str(minpc));
  set(handl(3,3),'UserData',1);
elseif strcmp(action,'coefs')
  minpc   = round(get(handl(3,6),'Value'));
  set(handl(2,4),'String','Now working on final PCR model')
  [t,p,b,ssq,eigs] = pcr(x,y,minpc);
  b = b((minpc-1)*ny+1:minpc*ny,:)';
  figure(handl(1,1))
  plot(b), hold on, plot(b,'o'), plot(zeros(nx,1),'-g'), hold off
  s = sprintf('Regression Coefficients in Final Model with %g pcs',minpc);
  title(s)
  xlabel('Variable Number')
  ylabel('Regression Coefficient')
  set(handl(3,3),'UserData',2);
  set(handl(4,12),'Visible','On');
  set(handl(2,4),'String','Click RESUME to Continue');
elseif strcmp(action,'qlims')
  if nargout > 5
    % Calculate qlim
    res = sum(((x - t*p').^2)');
    th1 = sum(eigs(minpc+1:min([mx nx])));
    th2 = sum((eigs(minpc+1:min([mx nx]))).^2);
    th3 = sum((eigs(minpc+1:min([mx nx]))).^3);
    h0 = 1 - ((2*th1*th3)/(3*th2^2));
    if h0 <= 0, h0 = 0.0001, end
    qlim = th1*(((1.65*sqrt(2*th2*h0^2)/th1) + 1 + th2*h0*(h0-1)/th1^2)^(1/h0));
    s1 = sprintf('The 95 Percent Q limit is %g',qlim);
    figure(handl(1,1))
    plot(1:mx,res,1:mx,res,'+',[1 mx],[qlim qlim],'--g')
    s = sprintf('Value of Q with 95 Percent Limit Based on %g PC Model',minpc);
    title(s)
    xlabel('Sample Number')
    ylabel('Value of Q')
  end
  figure(gui)
  set(handl(3,3),'UserData',3);
  set(handl(4,12),'Visible','On');
  set(handl(2,4),'String',[s1,' - Click RESUME to Continue']);
elseif strcmp(action,'tlims')
  if nargout > 6
    % Calculate t2lim
    if minpc > 1
      if mx > 300
        t2lim = (minpc*(mx-1)/(mx-minpc))*ftest(.05,minpc,300);
      else
        t2lim = (minpc*(mx-1)/(mx-minpc))*ftest(.05,minpc,mx-minpc);
      end
      s1 = sprintf('The 95 Percent T^2 limit is %g',t2lim);
      tvar = std(t);
	  tsqvals = sum((auto(t)').^2);
	  figure(handl(1,1))
	  plot(1:mx,tsqvals,1:mx,tsqvals,'+',[0 mx],[t2lim t2lim],'--g')
	  s = sprintf('Value of T^2 with 95 Percent Limit Based on %g PC Model',minpc);
      title(s)
      xlabel('Sample Number')
      ylabel('Value of T^2')
    else
	  tvar = std(t);
      s1 = 'T^2 not calculated when number of latent variables = 1';
      t2lim = 1.96^2;
    end
    figure(gui)
    set(handl(3,3),'UserData',4);
    set(handl(4,12),'Visible','On');
    set(handl(2,4),'String',[s1,' - Click RESUME to Continue']);
  end 
end