modlpred.m

PLS_Toolbox是用于故障检测与诊断方面的matlab工具箱

function [yprdn,resn,tsqn,scoresn] = modlpred(newx,modl,plots,q,w,lv,p)
%MODLPRED Predictions based on models created by MODLGUI
%  MODLPRED makes Y-block predictions based on an input X-block
%  and an existing regression model created using MODLGUI.
%  Inputs are the X-block data (newx) in the units of the
%  original data, the structure variable that contains the
%  regression model (modl), and an optional variable (plots)
%  which suppresses the plots when set to 0. Outputs are the
%  Y-block predictions (yprdn), residuals (resn), T^2 values
%  (tsqn), and the scores (scoresn).
%
%  NOTE: (newdata) will be scaled in MODLPRED using data
%  contained in (modl).
%
%I/O: [yprdn,resn,tsqn,scoresn] = modlpred(newx,modl,plots);
%
%MODLPRED can also make predictions based on an existing PLS model
%  constructed with the NIPALS algorithm from the PLS function.
%  Inputs are the matrix of predictor variables (newx),
%  the PLS model inner-relation coefficients (bin), the
%  x-block loadings (p), the y-block loadings (q), the
%  x-block weights (w), the number of  latent  variables to
%  use in prediction (lv), and an optional variable (plots)
%  which suppresses the plots when set to 0. Outputs are the
%  Y-block predictions (yprdn), residuals (resn), and the 
%  scores (scoresn). Note that T^2 are not calculated.
%
%I/O: [yprdn,resn,scoresn] = modlpred(newx,bin,p,q,w,lv,plots);
%
%See also: MODLGUI, PCAGUI, PCAPRO, PLS, XPLDST

%Copyright Eigenvector Research, Inc. 1998
%nbg, 11/98, 12/98

if nargin<2
  error('Error - insufficient number of inputs to MODLPRED')
end
if isa(modl,'struct')
  if nargin < 3
    plots = 1;
  end
  [mx,nx] = size(newx);
  if ~strcmp(class(modl),'struct')
    error('Input (modl) must be structure array from MODLGUI')
  end
  if nx~=size(modl.reg,1)
    error('X-block does not have correct number of variables')
  end
  switch modl.scale
  case 'mean'
    newx    = scale(newx,modl.meanx);
    yprdn   = rescale(newx*modl.reg,modl.meany);
  case 'auto'
    newx    = scale(newx,modl.meanx,modl.stdx);
    yprdn   = rescale(newx*modl.reg,modl.meany,modl.stdy);
  otherwise
    yprdn   = newx*modl.reg;
  end
  lvs     = size(modl.loads,2);

  if plots ~= 0
    plot(1:mx,yprdn,'+-')
    xlabel('Sample Number')
    ylabel('Predicted Value')
    title('New Sample Predictions') 
    pause
  end
  if nargout>1
    if strcmp(modl.name,'nip')
      scoresn = newx*modl.wts*inv(modl.loads'*modl.wts);
      resn    = (newx-scoresn*modl.loads').^2;  
    elseif strcmp(modl.name,'sim')
      scoresn = newx*modl.wts(:,1:lvs);
      for ii=1:size(scoresn,2)
        anorm = norm(modl.scores(:,ii));
        scoresn(:,ii) = scoresn(:,ii)*anorm;
      end
      resn    = (newx-scoresn*modl.loads').^2;
    elseif strcmp(modl.name,'pcr')
      scoresn = newx*modl.loads;
      resn    = (newx-scoresn*modl.loads').^2;
    else
      error('Model of unkown type')
    end
    if nx>1
      resn = sum(resn')';
    end
    ps     = 1./sqrt(diag(modl.scores'*modl.scores)/(modl.samps-1));
    tsqn   = (scoresn*diag(ps)).^2;
    tsqn   = sum(tsqn,2);
    if plots~=0
      plot(1:mx,resn,'+-r')
      hline(modl.reslim,'--g')
      xlabel('Sample Number')
      ylabel('Q Residual')
      title('Q Residuals with 95 Percent Limit')
      pause
      plot(1:mx,tsqn,'+-b')
      hline(modl.tsqlim,'--g')
      xlabel('Sample Number')
      ylabel('T^2')
      title('T^2 with 95 Percent Limit')
      pause
      plot(resn,tsqn,'+r')
      hline(modl.tsqlim,'--g'), vline(modl.reslim,'--g')
      s    = ' '; s = [s(ones(mx,1),:),int2str([1:mx]')];
      text(resn,tsqn,s)
      xlabel('Q Residual for New Sample')
      ylabel('T^2 for New Sample')
      title('T^2 vs Q for with 95 Percent Limits')
    end
  end
else
  if nargin < 6
    error('Error - insufficient number of inputs to MODLPRED')
  end
  if nargin < 7
    p     = 1;
  end
  if nargout > 3
    error('Error - T^2 not calculated for non-structure input')
  end
  bin     = modl; modl = plots; plots = p; p = modl; clear modl
  [mx,nx] = size(newx);
  [mq,nq] = size(q);
  [mw,nw] = size(w);
  that    = zeros(mx,lv);
  yprdn   = zeros(mx,mq);
  if lv>nw
    error(sprintf('Maximum number of latent variables exceeded (Max = %g)',nw));
  end
  x       = newx;
  for ii=1:lv
    that(:,ii) = x*w(:,ii);
    x  = x - that(:,ii)*p(:,ii)';
  end
  for ii=1:lv
    yprdn = yprdn + bin(1,ii)*that(:,ii)*q(:,ii)';
  end
  if plots ~= 0
    plot(1:mx,yprdn,'+-')
    xlabel('Sample Number')
    ylabel('Predicted Value')
    title('New Sample Predictions') 
    pause
  end
  if nargout>1
    tsqn = newx*w*inv(p'*w);  %tsqn = scoresn for NIPALS pieces input
    resn = (newx-tsqn*p').^2;  
    if nx>1
      resn = sum(resn')';
    end
    if plots~=0
      plot(1:mx,resn,'+-r')
      xlabel('Sample Number')
      ylabel('Q Residual')
      title('Q Residuals')
      pause
    end
  end
end