function [r2,mrse,termsselected,r2terms,mrseterms,coeff, r2adj]=mplregstep(X,Y,nterms)
%stepwise mpl regression. Orders regressors to variation explained
%and performs mpl regression accordingly
%Usage:
%[r2,mrse,termsselected,r2terms,mrseterms,coeff,r2adj]=mplregstep(X,Y,nterms)
%
%Inputs:
%X(n*m)         :m independent variables
%Y(n)           :dependend variable
%nterms         :number of terms to select for regression, can't be larger than m
%
%Outputs:
%r2             : variation explained
%mrse           : mean root squared error
%termsselected  : independent variables selected
%r2terms        : r2 of model after adding selected term
%mrseterms      : mrse of model after adding selected term
%coeff          : coefficients of constant + nterms in predicting model
%uses mplreg

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

n=size(X,1);
if nargin<3
    nterms=size(X,2);
else
    nterms=min(nterms,size(X,2));
end
termsselected=[];
for i=1:nterms
    [r2,mrse,best,rssbest]=mplregstep2(X,Y,termsselected);
    if length(r2)>1
        disp('Two regressors with equal r2, first one selected');
        r2terms(i)=r2(1);%if more than 1 with equal result, select first option
        r2adj(i)=1- (1-r2(1)).*( (n-1)/(n-nterms));
   
        mrseterms(i)=mrse(1);
        termsselected(i)=best(1);
    else
       r2terms(i)=r2;
        r2adj(i)=1- (1-r2).*( (n-1)/(n-nterms));
   
        mrseterms(i)=mrse;
        termsselected(i)=best;
    end
end
[r2,mrse,yf,res,coeff]=mplreg(X(:,termsselected),Y,'no');


%disp('selected columns in X: ');
%disp(terms)


function [r2,mrse,best,rssbest]=mplregstep2(X,Y,termsselected)

nrows=size(X,1);
N=size(X,2);
j=0;X1=[];X2=[];

%calculate residual sum of squares per column
if ~isempty(termsselected) 
    for i=1:N
        if ~isempty(find(i==termsselected))
            X1=[X1 X(:,i)];
        else
            j=j+1;
            X2=[X2 X(:,i)];
            cols(j)=i;
        end       
    end
else
    X2=X;
    cols=[1:1:N];
end

N2=size(X2,2);
for i=1:N2
    X3=[X1 X2(:,i)];
    [r2(i),mrse(i),yf,res,coeff]=mplreg(X3,Y,'no');
    rss(i)=sum(res.*res);
end
a=find(rss==min(rss));
best=cols(a);
rssbest=rss(a);
mrse=mrse(a);
r2=r2(a);

