plscvrg.m

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

function [press,cumpress,minlv,b,r,w,p,qlim,t2lim,tvar] = plscvrg(x,y,split,itern,lv,np,mc,gui,action)
%PLSCVRG SDEP PLS Cross Validation for use with MODLGUI

% Copyright
% Eigenvector Technologies
% 1995
global cumpress press minlv p q w t u bb ssqdif b r qlim t2lim itern 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 < 6
    np = 1;
  end
  if nargin < 7
    mc = 1;
  end
  press = zeros(split*itern*ny,lv);
  kk  = 0;
  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 k = 1:itern
    s = sprintf('Now working on iteration %g out of %g',k,itern);
	set(handl(2,4),'String',s)
    [z,inds] = sort(randn(mx,1));
    xx = x(inds,:);
    yy = y(inds,:);
    for i = 1:split
      kk = kk + 1;
      if mc ~= 0
        [calx,mnsx] = mncn([xx(1:ind(i,1)-1,:); xx(ind(i,2)+1:mx,:)]);
        testx = scale(xx(ind(i,1):ind(i,2),:),mnsx);
        [caly,mnsy] = mncn([yy(1:ind(i,1)-1,:); yy(ind(i,2)+1:mx,:)]);
        testy = scale(yy(ind(i,1):ind(i,2),:),mnsy);
      else
        calx = [xx(1:ind(i,1)-1,:); xx(ind(i,2)+1:mx,:)];
        testx = xx(ind(i,1):ind(i,2),:);
        caly = [yy(1:ind(i,1)-1,:); yy(ind(i,2)+1:mx,:)];
        testy = yy(ind(i,1):ind(i,2),:);
      end
      if ny > 1
        [p,q,w,t,u,b,ssqdif] = pls(calx,caly,lv);
        bbr = conpred(b,w,p,q,lv);
        for jj = 2:lv
	      i1 = (jj-1)*ny+1; i0 = i1-ny;
	      bbr(i1:jj*ny,:) = bbr(i1:jj*ny,:) + bbr(i0:(jj-1)*ny,:);
        end
      else
        bbr = simpls1(calx,caly,lv);
      end
      for j = 1:lv
        ypred = testx*bbr((j-1)*ny+1:j*ny,:)';
	    press((kk-1)*ny+1:kk*ny,j) = sum((ypred-testy).^2)';
      end
    end
	figure(handl(1,1))
    plot(sum(press(kk*ny-split*ny+1:kk*ny,:)))
    txt = sprintf('Cumulative PRESS for Iteration Number %g out of %g',k,itern);
    title(txt)
    xlabel('Number of Latent Variables')
    ylabel('PRESS')
    drawnow
  end  
  pause(2)
  cumpress = sum(press)/itern;
  figure(handl(1,1))
  plot(cumpress)
  title('Cumulative PRESS as a Function of Number of Latent Variables')
  xlabel('Number of Latent Variables')
  ylabel('PRESS')
  drawnow
  [a,minlv] = min(cumpress);
  figure(gui)
  s = sprintf('Minimum Cumulative PRESS is at %g LVs',minlv);
  if minlv>1
    s1 = ' Click slider to choose different number of LVs - then RESUME to continue';
    set(handl(3,6),'Max',minlv,'Value',minlv);
  else
    s1 = ' - Click RESUME to continue';
	set(handl(3,6),'Value',minlv);
	set(handl([3 5 6],6),'Visible','Off')
  end
  set(handl(4,6),'String',num2str(minlv))
  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(minlv));
  set(handl(3,3),'UserData',1);
elseif strcmp(action,'coefs')
  minlv    = round(get(handl(3,6),'Value'));
  set(handl(2,4),'String','Now working on final PLS model')
  [p,q,w,t,u,bb,ssqdif] = pls(x,y,minlv);
  b = conpred1(bb,w,p,q,minlv);
  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 LVs',minlv);
  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 > 4
    r = w*inv(p'*w)*inv(t'*t)*t';
  end
  if nargout > 7
    % Calculate qlim
    res = sum(((x - t*p').^2)');
    th1 = sum(res)/(mx - 1);
    th2 = (min([mx nx])-minlv)*((th1/(min([mx nx])-minlv))^2);
    th3 = (min([mx nx])-minlv)*((th1/(min([mx nx])-minlv))^3);
    h0 = 1 - ((2*th1*th3)/(3*th2^2));
    qlim = th1*(((2.33*sqrt(2*th2*h0^2)/th1) + 1 + th2*h0*(h0-1)/th1^2)^(1/h0));
    s1 = sprintf('The Approximate 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 Approximate 95 Percent Limit Based on %g LV Model',minlv);
    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 > 8
    % Calculate t2lim
    if minlv > 1
      if mx > 300
        t2lim = (minlv*(mx-1)/(mx-minlv))*ftest(.05,minlv,300);
      else
        t2lim = (minlv*(mx-1)/(mx-minlv))*ftest(.05,minlv,mx-minlv);
      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 LV Model',minlv);
      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