modlpred.m

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

function [ypred,res,tsq] = modlpred(newx,model,plots)
%MODLPRED Predictions based on models created by MODLMKER and MODLGUI
%  For regression models, this function makes new Y-block 
%  predictions using a new X-block and a compact model 
%  created by MODLMKER or MODLGUI. The outputs are the new vector 
%  (or matrix) of predicted variables (ypred), and for PLS/PCR
%  models the X-block residuals (res) and T^2 values (tsq). 
%
%  For PCA models, the function calculates the new scores (ypred),
%  residuals (res) and T^2 values (tsq).
%  
%  I/O format is  [ypred,res,tsq] = modlpred(newx,model,plots);
%
%  See also MODLMKER, MODLGUI, MODLRDER and MODLDEMO

%  Copyright
%  Barry M. Wise
%  Eigenvector Technologies
%  1994
%  Modified February 1995
%  Modified May 1995
[mx,nx] = size(newx);
[mm,nm] = size(model);
m = model;
if nargin < 3
  plots = 1;
end
if m(1,9) ~= 0
  if nx ~= m(1,8)
    error('New X block does not have correct number of variables')
  end
  if nargin == 2
    plots = 1;
  end
  if model(1,12) == 1
    me = 2+m(1,10)-1;
    b = model(2:me,1:m(1,8))';
    sx = newx;
    ypred = sx*b;
  elseif model(1,12) == 2
    me = 3+m(1,10)-1;
    b = model(3:me,1:m(1,8))';
    sx = scale(newx,m(2,1:m(1,8)));
    y = sx*b;
    ypred = rescale(y,m(2,m(1,8)+1:m(1,8)+m(1,10)));
  elseif model(1,12) == 3
    me = 4+m(1,10)-1;
    b = model(4:me,1:m(1,8))';
    sx = scale(newx,m(2,1:m(1,8)),m(3,1:m(1,8)));
    y = sx*b;
    my = m(2,m(1,8)+1:m(1,8)+m(1,10));
    sy = m(3,m(1,8)+1:m(1,8)+m(1,10));
    ypred = rescale(y,my,sy);
  else
    error('Model scaling not of known type')
  end
  if plots ~= 0
    plot(1:mx,ypred,'+',1:mx,ypred)
    xlabel('Sample Number')
    ylabel('Predicted Value')
    title('New Sample Predictions') 
    pause
  end
  if nargout > 1
    if m(1,11) ~= 3
      if m(1,11) == 1
        p = m(me+1:me+m(1,14),1:m(1,8));
        tvar = m(me+1:me+m(1,14),m(1,8)+1)';
	    t = sx*p';
	    tsq = sum(((scale(t,zeros(1,m(1,14)),tvar)).^2)');
	    res = sum(((sx - t*p)').^2);
      elseif m(1,11) == 2
        p = m(me+1:me+m(1,14),1:m(1,8));
        tvar = m(me+1:me+m(1,14),m(1,8)+1)';
	    w = m(me+1+m(1,14):me+2*m(1,14),1:m(1,8)); tsq = 0;
	    t = sx*w'*inv(p*w');
	    tsq = sum(((scale(t,zeros(1,m(1,14)),tvar)).^2)');
	    res = sum(((sx - t*p).^2)');
      end
	  if plots ~= 0
        plot(1:mx,res,1:mx,res,'+g',[0 mx],[m(1,15) m(1,15)],'--g')
        xlabel('Sample Number')
        ylabel('Q Residual')
        title('Q Residuals of New Samples with 95 Percent Limit')
        pause
        plot(1:mx,tsq,1:mx,tsq,'+g',[0 mx],[m(1,16) m(1,16)],'--g')
        xlabel('Sample Number')
        ylabel('T^2')
        title('T^2 for New Samples with 95 Percent Limit')
	    pause
	    plot(res,tsq,'+r'), hold on
	    z = axis;
	    axis(axis);
	    plot([z(1) z(2)],[m(1,16) m(1,16)],'--g')
	    plot([m(1,15) m(1,15)],[z(3) z(4)],'--g')
	    for i = 1:mx
	      text(res(i),tsq(i),int2str(i));
	    end
	    xlabel('Q Residual for New Sample')
	    ylabel('T^2 for New Sample')
	    title('T^2 versus Q Residual for New Samples with 95 Percent Limits')
	    hold off; axis('auto');
	  end
    else
      disp('Q and T^2 Plots not available for RR Models')
	  res = 0; tsq = 0;
    end
  end
else
  q = m(1,15); tsq = m(1,16); mo = m(1,7); ssq = zeros(m(1,14),4);
  if m(1,12) == 1
    loads = m(3:mm,1:m(1,8))';
    sx = newx;
	ssq(:,2) = m(2,1:m(1,14))'; 
  elseif model(1,12) == 2
    loads = m(4:mm,1:m(1,8))';
    sx = scale(newx,m(2,1:m(1,8)));
    ssq(:,2) = m(3,1:m(1,14))';
  elseif model(1,12) == 3
    loads = m(5:mm,1:m(1,8))';
    sx = scale(newx,m(2,1:m(1,8)),m(3,1:m(1,8)));
	ssq(:,2) = m(4,1:m(1,14))';
  end
  [ypred,res,tsq] = pcapro(sx,loads,ssq,q,tsq,plots,mo);
end