modlmkeg.m

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

function model = modlmkeg(x,y,gui,action)
%MODLMKEG Develops PCR, PLS and RR models.
% This  function is run by MODLGUI.
% It is used to develop PCR, PLS and ridge regression
% models with a variety of scaling and cross-validation
% options. The function is called through a graphical user
% interface. Inputs are the matrix of predictor variables (x)
% and predicted variable (y). The output is the regression
% model in compact form (model) that includes information
% about how the model was developed, scaling parameters, the final
% regression vector and a date and time stamp. The function
% also plots actual vs. fitted values and leverage vs.
% studentized residuals for outlier detection.
% I/O format is  model = modlmkeg(x,y,gui);
%
% See also MODLGUI, MODLPRED, MODLRDER and MODLDEMO

% Copyright Eigenvector Technologies 1995
global cumpress press minpc b rinv w p qlim t2lim
global split fys theta cvnum tvar handl
handl     = get(gui,'Userdata');
x         = get(handl(2,2),'Userdata');
y         = get(handl(2,3),'Userdata');

[smx,snx] = size(x);
[smy,sny] = size(y);

if get(handl(4,2),'Value')
  [ax,mx] = mncn(x); sx = ones(1,snx);
  [ay,my] = mncn(y); sy = ones(1,sny);
  sopt    = 2;
elseif get(handl(5,2),'Value')
  [ax,mx,sx] = auto(x);
  [ay,my,sy] = auto(y);
  sopt    = 3;
else
  ax = x; mx = 0; sx = 1;
  ay = y; my = 0; sy = 1;
  sopt    = 1;
end

if get(handl(3,3),'Value');
  ropt = 2;
elseif get(handl(4,3),'Value');
  ropt = 1;
else
  ropt = 3;
end
if ropt~=3
  if get(handl(3,5),'Value');
    copt = 1;
  elseif get(handl(4,5),'Value');
    copt = 2;
  elseif get(handl(5,5),'Value');
    copt = 3;
  else
    copt = 4;
  end
else
  if get(handl(9,5),'Value');
    copt = 1;
  else get(handl(8,5),'Value');
    copt = 2;
  end
end

if strcmp(action,'final')~=1
  if ropt == 1
    pc = round(get(handl(3,6),'Value'));
    if copt == 1
	  split = round(get(handl(3,7),'Value'));
      [press,cumpress,minpc,b,p,qlim,t2lim,tvar] = pcrcvg(ax,ay,split,pc,1,0,gui,action);
    elseif copt == 2
      [press,cumpress,minpc,b,p,qlim,t2lim,tvar] = pcrcv1g(ax,ay,pc,1,0,gui,action);
    elseif copt == 3
	  split = round(get(handl(3,7),'Value'));
      [press,cumpress,minpc,b,p,qlim,t2lim,tvar] = pcrcvbg(ax,ay,split,pc,1,0,gui,action); 
    elseif copt == 4
	  split = round(get(handl(3,7),'Value'));
      cvnum = round(get(handl(3,8),'Value'));
      [press,cumpress,minpc,b,p,qlim,t2lim,tvar] = pcrcvrg(ax,ay,split,cvnum,pc,1,0,gui,action);
    end
  elseif ropt == 2
    pc = round(get(handl(3,6),'Value'));
    if copt == 1
	  split = round(get(handl(3,7),'Value'));
      [press,cumpress,minpc,b,rinv,w,p,qlim,t2lim,tvar] = plscvg(ax,ay,split,pc,1,0,gui,action);
    elseif copt == 2
      [press,cumpress,minpc,b,rinv,w,p,qlim,t2lim,tvar] = plscv1g(ax,ay,pc,1,0,gui,action);
    elseif copt == 3
	  split = round(get(handl(3,7),'Value'));
      [press,cumpress,minpc,b,rinv,w,p,qlim,t2lim,tvar] = plscvbg(ax,ay,split,pc,1,0,gui,action); 
    elseif copt == 4
	  split = round(get(handl(3,7),'Value'));
      cvnum = round(get(handl(3,8),'Value'));
      [press,cumpress,minpc,b,rinv,w,p,qlim,t2lim,tvar] = plscvrg(ax,ay,split,cvnum,pc,1,0,gui,action);
    end
  else
    tmax = round(get(handl(3,9),'Value')*100)/100;
    divs = get(handl(3,10),'Value');
    if copt == 1
      [b,theta] = ridgeg(ax,ay,tmax,divs,0,gui,action);
    elseif copt == 2
	  split = round(get(handl(3,11),'Value'));
      [b,theta,cumpress] = ridgecg(ax,ay,tmax,divs,split,gui,action);
    end
  end
elseif strcmp(action,'final')
  minpc = round(get(handl(3,6),'Value'));
  if get(handl(3,3),'UserData') ==4
    if ropt ~= 3
      if sny == 1
        if sopt ~= 3
          if copt == 4
	        rmsecv = sqrt(cumpress(minpc)/(split*cvnum));
	      else
            rmsecv = sqrt(cumpress(minpc)/smx);
	      end
        else
          if copt == 4
	        rmsecv = sqrt(cumpress(minpc)*(sy.^2)/(split*cvnum));
	      else
            rmsecv = sqrt(cumpress(minpc)*(sy.^2)/smx);
	      end
	    end	
      else
	    rmsecv = zeros(sny,1);
	    [mpress,npress] = size(press);
        for yvar = 1:sny
	      rmsecv(yvar) = sum(press(yvar:mpress,minpc));
  	    end	
	    if sopt ~= 3
          if copt == 4
	        rmsecv = sqrt(rmsecv/(split*cvnum));
	      else
            rmsecv = sqrt(rmsecv/smx);
	      end
        else
          if copt == 4
	        rmsecv = sqrt(rmsecv.*(sy.^2)'/(split*cvnum));
	      else
           rmsecv = sqrt(rmsecv.*(sy.^2)'/smx);
	      end
		end
	  end	
    end
    rmsf = 0;
    if ropt == 3
      if copt == 1
        rmsf = 1;
      else
	    rmsecv = sqrt(min(cumpress)*(sy^2)/smx);
      end
    end
    if rmsf == 0
	  figure(gui);
      if sny == 1
        s1 = sprintf('Root Mean Square Error of Cross Validation = %g',rmsecv);
      else
        s1 = ['Root Mean Square Errors of Cross Validation'];
	    s2 = sprintf('for the %g Y-block Variables are:',sny);
	    s3 = num2str(rmsecv);
		s1 = [s1,s2,s3];
      end
    end
    % Plot fitted versus true value
    figure(handl(1,1));
    fy    = ax*b;
    fys   = rescale(fy,my,sy);
    plot(y,fys,'+')
    for i = 1:smx
      for j = 1:sny
        s   = int2str(i);
        text(y(i,j),fys(i,j),s)
      end
    end
    z = axis;
    hold on, plot(z(1:2),z(1:2),'-g'), hold off
    xlabel('Actual Value')
    ylabel('Predicted Value'), drawnow
    set(handl(3,3),'UserData',5);
    set(handl(4,12),'Visible','On');
    set(handl(2,4),'String',[s1,' - Click RESUME to Continue']);	
  elseif get(handl(3,3),'UserData') == 5
    % Calculate RMSEC
    rmsec  = sqrt(sum((fys-y).^2)/smx);
    if sny == 1
      s = sprintf('Root Mean Square Error of Calibration is %g',rmsec);
      set(handl(2,4),'String',s)
    else
      s  = sprintf('for the %g Y-block Variables are: ',sny);
      s2 = sprintf('Root Mean Square Errors of Calibration ');
      s3 = num2str(rmsec');
      set(handl(2,4),'String',[s2, s, s3])
    end
    set(handl(3,3),'UserData',6);
    set(handl(4,12),'Visible','On');
	s = get(handl(2,4),'String');
    set(handl(2,4),'String',[s,' - Click RESUME to Continue']);	
  elseif get(handl(3,3),'UserData') == 6
    % Calculate R+
    if ropt == 1
      t    = ax*p;
      rinv = p(:,1:minpc)*inv(t(:,1:minpc)'*t(:,1:minpc))*t(:,1:minpc)';
    elseif ropt == 2
    % rinv calculated in pls routine
    else
      rinv = inv(ax'*ax + eye(snx)*theta)*ax';
    end
    lev   = zeros(smx,1);
    for i = 1:smx
      lev(i) = ax(i,:)*rinv(:,i);
    end
    res   = y - fys;
    sigcv = zeros(1,sny);
    tres  = zeros(smx,sny);
    ecv   = zeros(smx,sny);
    for i = 1:sny
      sigcv(i)  = sqrt(inv(smx-1)*res(:,i)'*(((ones(smx,1)-lev).^(-2)).*res(:,i)));
      tres(:,i) = res(:,i)./(sigcv(i)*sqrt(ones(smx,1)-lev));
      ecv(:,i)  = res(:,i).*(ones(smx,1)./(ones(smx,1)-lev));
    end
    rmsecv = sqrt(sum(ecv.^2)/smx);
    if sny == 1
      s1 = ['Root Mean Square Error of Leave One Out '];
      s  = sprintf('Cross Validation is %g',rmsecv);
      s1 = [s1,s];
    else
      s1 = ['Root Mean Square Errors of Leave One Out '];
      s  = sprintf('Cross Validation for the %g Y-block variables are:',sny);
      s2 = num2str(rmsecv');
      s1 = [s1, s, s2];
    end
	figure(handl(1,1))
    plot(lev,tres,'+')
    for i = 1:smx
      for j = 1:sny
        s = [' ',int2str(i)]; 
        text(lev(i),tres(i,j),s);
      end
    end
    xt = get(gca,'XTick');
    [mxt,nxt] = size(xt);
    hold on, plot([xt(1) xt(nxt)],[0 0],'-g'), hold off
    xlabel('Leverage')
    ylabel('Studentized Residuals')

    if snx+sny > 16
      model = zeros(1,snx+sny); sm = snx+sny;
    else
      model = zeros(1,16); sm = 16;
    end
   %if nargin > 1
      model(1,1:13) = [clock smx snx smy sny ropt sopt copt];
      if ropt == 3
        model(1,14) = theta;
      else
        model(1,14:16) = [minpc qlim t2lim];
      end
    %else
    %  model(1,1:14) = [clock smx snx smy q tsq sopt 0 spca];
    %end
    if sopt ~= 1
      model                  = [model; zeros(1,sm)];
      model(2,1:snx)         = mx;
      model(2,snx+1:snx+sny) = my;
      if sopt == 3
        model                  = [model; zeros(1,sm)];
        model(3,1:snx)         = sx;
        model(3,snx+1:snx+sny) = sy;
      end
    end
    reg = zeros(sny,sm);
    reg(1:sny,1:snx) = b';
    model = [model; reg];
    if ropt == 1
	  model = [model; [p' tvar' zeros(minpc,sm-snx-1)]];
    end
    if ropt == 2
  	  model = [model; [p' tvar' zeros(minpc,sm-snx-1)]; [w' zeros(minpc,sm-snx)]];
    end
    set(handl(3,12),'Userdata',model);
    s1  = [s1,' - Model Complete - Click QUIT to Save MODEL and Exit'];
    set(handl(2,4),'String',s1);
  end
end %end of elseif action == 'final'