plscvg.m

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

function [press,cumpress,minlv,b,r,w,p,qlim,t2lim,tvar] = plscvg(x,y,split,lv,np,mc,gui,action)
%PLSCVG Cross validation for PLS models used with MODLGUI

% Copyright
% Eigenvector Technologies
% 1995

global cumpress press minlv p q w t u bb ssqdif b r qlim t2lim 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,lv);
  for i = 1:split
    s = sprintf('Now working on test set %g out of %g',i,split);
    set(handl(2,4),'String',s);
    ind = zeros(mx,1);
    count = 0;
    for j = 1:mx
      test = round((j+i-1)/split) - ((j+i-1)/split);
      if test == 0
        ind(j,1) = 1;
        count = count + 1;
      end
    end
    [a,b] = sort(ind);
    newx = x(b,:);
    newy = y(b,:);
    if mc ~= 0
      [calx,mnsx] = mncn(newx(1:mx-count,:));
      testx = scale(newx(mx-count+1:mx,:),mnsx);
      [caly,mnsy] = mncn(newy(1:mx-count,:));
      testy = scale(newy(mx-count+1:mx,:),mnsy);
    else
      calx = newx(1:mx-count,:);
      testx = newx(mx-count+1:mx,:);
      caly = newy(1:mx-count,:);
      testy = newy(mx-count+1:mx,:);
    end
    if ny > 1
      [p,q,w,t,u,b,ssqdif] = pls(calx,caly,lv);
      bbr = conpred(b,w,p,q,lv);
      for k = 2:lv
	    i1 = (k-1)*ny+1; i0 = i1-ny;
	    bbr(i1:k*ny,:) = bbr(i1:k*ny,:) + bbr(i0:(k-1)*ny,:);
      end
    else
      bbr = simpls1(calx,caly,lv);
    end
    for j = 1:lv
      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 Latent Variables')
    ylabel('PRESS')
    drawnow
  end
  pause(2)
  cumpress = sum(press);
  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