function [r2,rmse,yf,res,coeff, re, r2adj,cimin,cimax,index_outliers,index_not_outliers]=mplreg(X,Y,display, order, intercept, p)
%Multiple linear regression
%
%USAGE:
%[r2,rmse,yf,res,coeff, re, r2adj]=mplreg(X,Y,display, order, intercept)
%
%INPUT:
%X [n*m]            : independent variables
%Y [n]              : dependent variable
%display            : 'yes' for display
%intercept          : 1/0
%p                  : p value for confidence interval. Default = 0.9
%
%OUTPUT:
%r2                 : explained variation
%rmse               : root mean squared error
%yf    [n]          : fitted values
%res   [n]          : residuals
%coeff [m + 1]      : regression coefficients
%re                 : relative error
%r2adj              : adjusted r2
%cimin              : lower confidence interval
%cimax              : upper confidence interval
%index_outliers     : index of elements that are outside confidence interval
%index_not_outliers : index of elements that are within confidence interval
%see also mplpoly for polynomal fit for 1 X variate
%
%Tom Schut, March 2005

if nargin <3
    display='yes';
end
if nargin <4
    order=0;
end
if nargin <5
    intercept=1;
end
if nargin <6
    p=0.9;
end

a=find(~isnan(mean(X,2)) & ~isnan(Y));
b=find(isnan(mean(X,2)) | isnan(Y));
X=X(a,:); 
Y=Y(a);

Xord=[];
for ord=2:order
    Xord=[Xord X.^ord];
end
X=[X Xord];



nd=size(X,1);
N=size(X,2);
   

if nd==0 | length(Y)==0
    r2=NaN;
    r2adj=NaN;
    rmse=NaN;
    yf=NaN;
    res=NaN;
    coeff=[NaN, NaN];
    re=NaN;
%elseif(N-nd)<0    
    whos X
    whos Y
    error('nr of regressors equal or larger than nr of observations');
elseif(size(X,1) ~= size(Y,1))
    error('X and Y should have equal nr of rows');
%elseif length(b)>0
%    disp(sprintf('a total of %d records with NaN found and ignored',length(b)));    
end

y=Y;
if intercept
    x=[ones(nd,1),X];
else
    x=X;
end    
coeff=x\y;
yf=x*coeff;

se=(y-yf).^2;
res=(y-yf);
rmse=sqrt( mean(res.*res) );
re=rmse ./ mean(y);
r2=1-sum(se)/sum((y-mean(y)).^2);
r2adj=1- (1-r2).*( (nd-1)/(nd-N));

if strcmp(display,'yes')
    disp(sprintf('R2 = %6.4f, R2adj = %6.4f, rmse = %6.2f, re= %6.2f, N= %6.0f',r2,r2adj, rmse,re,nd));
    if intercept
        s=sprintf('Model fitted: Y= %5.6f',coeff(1));
        for i=1:N
            s=strcat(s,sprintf(' + %5.6fX%1.0f',coeff(i+1),i));
        end
    else
        s=sprintf('Model fitted: Y= %5.6fX',coeff(1));
    end
    disp(s);
end

%find confidence interval of new observation
%--------------------------------------------------------------------
    SSX=(y-se).^2;
    MSE=sum(se)./nd;
    S=sqrt(MSE.*(1/nd + 1/1 + SSX./sum(SSX)));
    
    t=tinv(0.5+p/2,nd-2);
    
    %the fitted line is the 1:1 line, where x=y
    cimin(a,1)=y-t.*S;
    cimax(a,1)=y+t.*S;
%--------------------------------------------------------------------

outlier_value(a,1)=abs(res)-t.*S;
outlier_value(b,1)=NaN;

YF(b,1)=NaN;
YF(a,1)=yf;
yf=YF;
RES(b,1)=NaN;
RES(a,1)=res;
res=RES;
cimax(b,1)=NaN;
cimin(b,1)=NaN;    

index_outliers = find(outlier_value > 0);
index_not_outliers = find(outlier_value <= 0);