function result=predict(x,y,inmodel)
%PREDICT is a function which computes regression results for new data based on
% the output from a RPCR or RSIMPLS analysis.
%
% I/O: result=predict(x,y,inmodel)
%
% Required input arguments:
% x: Data matrix of the explanatory variables of the new data
% y: Data matrix of the response variables of the new data
% inmodel: output from a RPCR or RSIMPLS model
%
% The output of PREDICT is a structure containing:
%
% result.fitted : Fitted response values of the new data
% result.res : Residuals of the new data
% result.sd : Score distances of the new data
% result.od : Orthogonal distances of the new data
% result.resd : Residual distances of the new data
% result.cutoff : Cutoff values for the score (result.cutoff.sd), orthogonal
% (result.cutoff.od) and residual distances (result.cutoff.resd).
% We use 0.99 quantiles of the chi-squared distribution.
% result.flag : The observations whose score distance is larger than
% 'result.cutoff.sd' receive a flag 'result.flag.sd' equal
% to zero (good leverage points). Otherwise 'result.flag.sd'
% is equal to one.
% The components 'result.flag.od' and 'result.flag.resd' are
% defined analogously, and determine the orthogonal outliers,
% resp. the bad leverage points/vertical outliers.
% The observations with 'result.flag.od' and 'result.flag.resd'
% equal to zero, can be considered as calibration outliers and receive
% 'result.flag.all' equal to zero. The regular observations and the good leverage
% points have 'result.flag.all' equal to one.
% result.rmsep : the root mean squared error of the non-outlying data
% result.class : the name of the method used: 'RPCR' or 'RSIMPLS'
% result.classic : if there was a classical output from rpcr or rsimpls,
% all computations based on these classical results are
% done as well. Note that the root mean squared error is
% then based on the non-outlying data from the robust
% analysis, in order to compare the rmsep values on the
% same data set.
%
% Example:
% result.pcr = rpcr(Xtrain,ytrain,'k',3,'plots',0);
% result = predict(Xtest,ytest,result.pcr);
%
% Written by S. Verboven, M. Hubert
% Last Revision: 20/05/2008
X=x;
Y=y;
[n,p]=size(X);
[n,q]=size(Y);
%Fitted values and residuals of the new data
result.fitted=X*inmodel.slope + repmat(inmodel.int,n,1);
result.res=Y-result.fitted;
%Computing distances
if strcmp(inmodel.class,'RPCR')
XRc=X-repmat(inmodel.robpca.M,n,1);
result.T=XRc*inmodel.robpca.P;
result.sd=sqrt(libra_mahalanobis(result.T,zeros(size(result.T,2),1),'cov',inmodel.robpca.L))';
Xtilde=result.T*inmodel.robpca.P';
Rdiff=XRc-Xtilde;
for i=1:n
result.od(i,1)=norm(Rdiff(i,:));
end
if q >1
cen=zeros(q,1)';
result.resd=sqrt(libra_mahalanobis(result.res,cen,'cov',inmodel.mcdreg.cov))';
else
result.resd=result.res/inmodel.lts.scale;
end
%cutoffs
quan=chi2inv(0.99,inmodel.k);
result.cutoff.sd=quan^0.5;
if inmodel.k~=rank(X)
[m,s]=unimcd(inmodel.od.^(2/3),inmodel.h);
result.cutoff.od = sqrt(norminv(0.99,m,s).^3);
else
result.cutoff.od=0;
end
result.cutoff.resd=sqrt(chi2inv(0.99,q));
elseif strcmp(inmodel.class,'RSIMPLS')
XRc=X-repmat(inmodel.robpca.M(1:p),n,1);
result.T=XRc*inmodel.weights.r;
result.sd=sqrt(libra_mahalanobis(result.T,inmodel.Tcenter,'cov',inmodel.Tcov))';
Xtilde=result.T*inmodel.weights.p';
Rdiff=XRc-Xtilde;
for i=1:n
result.od(i,1)=norm(Rdiff(i,:));
end
if q >1
cen=zeros(q,1)';
result.resd=sqrt(libra_mahalanobis(result.res,cen,'cov',inmodel.cov))';
else
result.resd=result.res/sqrt(inmodel.cov);
end
%cutoffs
quan=chi2inv(0.99,inmodel.k);
result.cutoff.sd=quan^0.5;
if inmodel.k~=rank(X)
[m,s]=unimcd(inmodel.od.^(2/3),inmodel.h);
result.cutoff.od = sqrt(norminv(0.99,m,s).^3);
else
result.cutoff.od=0;
end
result.cutoff.resd=sqrt(chi2inv(0.99,q));
end
% Defining flags
result.flag.od=(result.od<=result.cutoff.od);
result.flag.sd=(result.sd<=result.cutoff.sd);
result.flag.resd=(abs(result.resd)<=result.cutoff.resd);
result.flag.all=result.flag.od & result.flag.resd;
% Rmsep based on data with 'result.flag.resd = 1'
N=sum(result.flag.resd);
if q>1
result.rmsep=sqrt(1/(N*q)*sum(sum(result.res(result.flag.resd==1).^2,2)));
else
result.rmsep=sqrt(1/N*sum((result.res(result.flag.resd==1)).^2));
end
result.class=inmodel.class;
%In case the classical output is also given
if isstruct(inmodel.classic) && strcmp(inmodel.class,'RPCR')
% fitted values, residuals, distances
result.classic.fitted=X*inmodel.classic.slope + repmat(inmodel.classic.int,n,1);
result.classic.res=Y-result.classic.fitted;
XRc=X-repmat(inmodel.classic.cpca.M,n,1);
result.classic.T=XRc*inmodel.classic.cpca.P;
result.classic.sd=sqrt(libra_mahalanobis(result.classic.T,zeros(size(result.classic.T,2),1),'cov',inmodel.classic.cpca.L))';
Xtilde=result.classic.T*inmodel.classic.cpca.P';
Rdiff=XRc-Xtilde;
for i=1:n
result.classic.od(i,1)=norm(Rdiff(i,:));
end
if q >1
cen=zeros(q,1)';
result.classic.resd=sqrt(libra_mahalanobis(result.classic.res,cen,'cov',inmodel.classic.cov))';
else
result.classic.resd=result.classic.res/sqrt(inmodel.classic.cov);
end
% cutoffs and flags
quan=chi2inv(0.99,inmodel.k);
result.classic.cutoff.sd=quan^0.5;
if inmodel.k~=rank(X)
[m,s]=unimcd(inmodel.classic.od.^(2/3),inmodel.h);
result.classic.cutoff.od = sqrt(norminv(0.99,m,s).^3);
else
result.classic.cutoff.od=0;
end
result.classic.cutoff.resd=sqrt(chi2inv(0.99,q));
result.classic.flag.od=(result.classic.od<=result.classic.cutoff.od);
result.classic.flag.sd=(result.classic.sd<=result.classic.cutoff.sd);
result.classic.flag.resd=(abs(result.classic.resd)<=result.classic.cutoff.resd);
result.classic.flag.all=result.classic.flag.od & result.classic.flag.resd;
% classical rmsep of the good observations (from the robust analysis)
N=sum(result.flag.resd);
if q>1
result.classic.rmsep=sqrt(1/(N*q)*sum(sum(result.classic.res(result.flag.resd==1).^2,2)));
else
result.classic.rmsep=sqrt(1/N*sum((result.classic.res(result.flag.resd==1)).^2));
end
result.classic.class=inmodel.classic.class;
elseif isstruct(inmodel.classic) && strcmp(inmodel.class,'RSIMPLS')
% fitted values, residuals, distances
result.classic.fitted=X*inmodel.classic.slope + repmat(inmodel.classic.int,n,1);
result.classic.res=Y-result.classic.fitted;
XRc=X-repmat(inmodel.classic.M(1:p),size(X,1),1);
result.classic.T=XRc*inmodel.classic.weights.r;
result.classic.sd=sqrt(libra_mahalanobis(result.classic.T,zeros(size(result.classic.T,2),1),'cov',inmodel.classic.Tcov))';
Xtilde=result.classic.T*inmodel.classic.weights.p';
Rdiff=XRc-Xtilde;
for i=1:n
result.classic.od(i,1)=norm(Rdiff(i,:));
end
if q >1
cen=zeros(q,1)';
result.classic.resd=sqrt(libra_mahalanobis(result.classic.res,cen,'cov',inmodel.classic.cov))';
else
result.classic.resd=result.classic.res/sqrt(inmodel.classic.cov);
end
% cutoffs and flags
quan=chi2inv(0.99,inmodel.k);
result.classic.cutoff.sd=quan^0.5;
if inmodel.k~=rank(X)
[m,s]=unimcd(inmodel.classic.od.^(2/3),inmodel.h);
result.classic.cutoff.od = sqrt(norminv(0.99,m,s).^3);
else
result.classic.cutoff.od=0;
end
result.classic.cutoff.resd=sqrt(chi2inv(0.99,q));
result.classic.flag.od=(result.classic.od<=result.classic.cutoff.od);
result.classic.flag.sd=(result.classic.sd<=result.classic.cutoff.sd);
result.classic.flag.resd=(abs(result.classic.resd)<=result.classic.cutoff.resd);
result.classic.flag.all=result.classic.flag.od & result.classic.flag.resd;
% classical rmsep of the good observations (from the robust analysis)
N=sum(result.flag.resd);
if q>1
result.classic.rmsep=sqrt(1/(N*q)*sum(sum(result.classic.res(result.flag.resd==1).^2,2)));
else
result.classic.rmsep=sqrt(1/N*sum((result.classic.res(result.flag.resd==1)).^2));
end
result.classic.class=inmodel.classic.class;
end