
function [S,b,r,pval]=regresion(X,Y,method,plotflag)
% [S,b,r,pval,MSE,N,SEm]=regresion(X,Y,method) % have commented out MSE and
% beyod
% Computes a Regression Line that intercepts zero
% label should be a string with the name of the variable
% gives the slope and the Mean-squared error
% Athird entry is needed to select the method:
% 0= Least Squares solution, No intercep (forced to cross (0,0))
% 1= Least Squares solution
% 2= Both
% plotflag==1; make plot, 0 do not show plot
% from http://www.mathworks.com/matlabcentral/newsreader/view_thread/158607
% Originally created under /Users/Mau/Documents/FALL2010/DataAnalysis_Gille/Project3/QscatvsWindsat

fts=14; %fontsize for legend

% Ignore -Inf (aded agust 19)
X(isinf(X))=NaN;
Y(isinf(Y))=NaN;

% Ignore NaNs (added Jan 03 2011)
ny=isnan(Y);
nx=isnan(X);

X(ny)=NaN;
Y(nx)=NaN;

x=X(~isnan(X));
y=Y(~isnan(Y));

[r, pval]=nancorr(x,y);%,'type','pearson','rows','complete');

% Compute number of obs and standard error for the slope (Feb 9, 2015)
% and do not compute of N=1
[nrows, ncolumns]=size(unique([x y],'rows'));
N=nrows;

if N<2
    N=N;
    S=NaN;
    b=NaN;
    r=NaN;
    pval=NaN;
    MSE=NaN;
    SEm=NaN;
elseif N>=2
    
    if method==0
        % No intercept Linear Model (Least-squares solution)
        % S = x\y; % slope, x = A\b returns a least squares solution for Ax=b, when a is not squared))
        [S,se_S,MSE] = lscov(x,y); % lscov allows to compute the mean-squared error
        b=0;
        
        yhat = x*S; % model
        EVR = var(yhat)/var(y); % "Explained" Variance Ratio
        plot(x,yhat,'k');
%         e = y-yhat; % error
%         MSE = mse(e); % mean-squared error
        
        % Standard error for the slope:
        % If you need to calculate the standard error of the slope (SE) by hand, use the following formula:
        % SE = sb1 = sqrt [ ?(yi - ?i)2 / (n - 2) ] / sqrt [ ?(xi - x)2 ]
        
        res = y - yhat;
        resx=x-mean(x);
        SEm=sqrt(sum(res.^2)/(N-2))./sqrt(sum(resx.^2));
        
        mRMSE=nanmean(sqrt(nanmean(res.^2)));
        
        plot(x,y,'o','color','k','markersize',6,'markerfacecolor',[0.8 0.8 0.8]);
        hold on
        plot(x,yhat,'k','linewidth',2);
        legend(['r = ', num2str(r, '%6.2f'),', pval = ', num2str(pval, '%6.2f')],...
            [' y =',num2str(S, '%6.4f'),' x, RMSE = ',num2str(mRMSE,'%6.4f' ),... MSE = ',num2str(MSE,'%6.4f' ),...
            ', N = ',num2str(N,'%6.0f' )],'location','S');
        grid on
        %hold off
        
    elseif method==1
        
        % Calculate fit parameters
        [p,ErrorEst] = polyfit(x,y,1);
        S=p(1); % Slope;
        b=p(2); % Ordenada al origen
        
        % Evaluate the fit
        y_fit = polyval(p,x,ErrorEst);
        % % Plot the data and the fit
        % plot(x,y_fit,'-',x,y,'+');
        
        EVR = var(y_fit)/var(y); % "Explained" Variance Ratio
        
        % % Annotate the plot
        % legend('Polynomial Model','Data','Location','NorthWest');
        % xlabel('Census Year');
        % ylabel('Population (millions)');
        
        % Calculate residuals
        res = y - y_fit;
%        % MSE = sqrt(sum(res.^2))./length(x(~isnan(x)));
%         MSE = mse(res);
        
        % Standard error for the slope:
        % If you need to calculate the standard error of the slope (SE) by hand, use the following formula:
        % SE = sb1 = sqrt [ ?(yi - ?i)2 / (n - 2) ] / sqrt [ ?(xi - x)2 ]
        
        resx=x-mean(x);
        SEm=sqrt(sum(res.^2)/(N-2))./sqrt(sum(resx.^2));
        
        % figure, plot(x,res,'+')
        % title('Residuals for the Linear Model')
        
        if plotflag==1
        % plot(x,y,'o','markersize',4,'markerfacecolor','b');
        plot(x,y,'o','color','k','markersize',6,'markerfacecolor','none'); %[0.3 0.3 0.3]);
        hold on
        plot(x,y_fit,'color',[0.2 0.2 0.2],'linewidth',2);
        
%         legend(['Correlation: r = ', num2str(r, '%6.2f'),', pval = ', ...
%             num2str(pval, '%6.2f'),', N=',num2str(N),' '],'location','NE');
        
%         legend(['r = ', num2str(r, '%6.2f'),', pval = ', ...
%             num2str(pval, '%6.3f')],[' y =',num2str(p(1), '%6.3f'),' x + ',... %num2str(p(1), '%6.5f')
%             num2str(p(2), '%6.3f')], ... %', MSE = ',num2str(MSE,'%6.2f' ), ', N=',num2str(N),' '],
%             'location','best');
        
        % legend in scientific notation
        legend(['r = ', num2str(r, '%6.2f'),', pval = ', ...
            num2str(pval, '%6.3f')],[' y =',num2str(p(1), '%6.1e'),' x + ',... %num2str(p(1), '%6.5f')
            num2str(p(2), '%6.1e')], ... %', MSE = ',num2str(MSE,'%6.2f' ), ', N=',num2str(N),' '],
            'location','SE'); %'best');
        
        grid on
        % hold off
        
        end
        
        % %All intercept Constant Model
        % y00 = mean(y); % 0.022617
        % e00 = y - mean(y) ;% error
        % MSE00 = mse(e00); % 8.0045
        %
        % R2 = 1 - MSE/MSE00; % R^2 = 0.76222 (So-called EVR)
        
    elseif method==2
        % No intercept linear model (least squares solution)
        [So,se_S,MSE] = lscov(x,y); % lscov allows to compute the mean squared error
        
        yhat = x*So; % model
        EVR = var(yhat)/var(y); % "Explained" Variance Ratio
        
        % % plot(x,y,'o','markersize',4,'markerfacecolor','b');
        plot(x,y,'o','color',[0.1 0.1 0.1],'markersize',7.5); %,'markerfacecolor','transparent');% [0.1 0.1 0.1]);
        hold on
        % plot(x,yhat,'r','linewidth',2);
        plot(x,yhat,'k','linewidth',2);
        hold on
        grid on
        
        % Linear regression (least squares solution)
        % Calculate fit parameters
        [p,ErrorEst] = polyfit(x,y,1);
        S=p(1);
        b=p(2);
        % Evaluate the fit
        y_fit = polyval(p,x,ErrorEst);
        
        % Calculate residuals
        res = y - y_fit;
%         % MSE = sqrt(sum(res.^2))./length(x(~isnan(x)));
%         MSE1 = mse(res);      
        
        % Standard error for the slope:
        % If you need to calculate the standard error of the slope (SE) by hand, use the following formula:
        % SE = sb1 = sqrt [ ?(yi - ?i)2 / (n - 2) ] / sqrt [ ?(xi - x)2 ]
        
        resx=x-mean(x);
        SEm=sqrt(sum(res.^2)/(N-2))./sqrt(sum(resx.^2));
        
        
        % plot(x,y_fit,'color',[0.2 0.8 0.2],'linewidth',2);
        plot(x,y_fit,'color',[0.4 0.4 0.4],'linewidth',2);
        h=legend(['Correlation: r = ', num2str(r, '%6.2f'),...
            ', pval = ', num2str(pval, '%6.3f')],[' y =',num2str(So, '%6.3f'),...
            ' x, MSE = ',num2str(MSE,'%6.2f' )],[' y =',num2str(p(1), '%6.3f'),...
            ' x + ',num2str(p(2), '%6.3f')],... % ' , MSE = ',num2str(MSE1,'%6.2f' )],...
            'location','NW');
        set(h,'fontsize',fts);
        
        % legend boxoff 
        hold off
        
        % %All intercept Constant Model
        % y00 = mean(y); % 0.022617
        % e00 = y - mean(y) ;% error
        % MSE00 = mse(e00); % 8.0045
        %
        % R2 = 1 - MSE/MSE00; % R^2 = 0.76222 (So-called EVR)
        
    elseif method==3
        % No intercept Linear Model (Least-squares solution)
        % S = x\y; % slope, x = A\b returns a least squares solution for Ax=b, when a is not squared))
        % [S,se_S,MSE] = lscov(x,y); % lscov allows to compute the mean-squared error
        
        % yhat = x*S; % model
        % EVR = var(yhat)/var(y); % "Explained" Variance Ratio
        % plot(x,yhat,'k');
        % e = y-yhat; % error
        % MSE = mse(e); % mean-squared error
        
        plot(x,y,'o','color','k','markersize',4,'markerfacecolor',[0.8 0.8 0.8]);
        hold on
        % plot(x,yhat,'k','linewidth',2);
        legend(['Correlation: r = ', num2str(r, '%6.2f'),', pval = ', num2str(pval, '%6.2f')],...
            [' y =',num2str(S, '%6.2f'),' x,  MSE = ',num2str(MSE,'%6.2f' )],'location','NE');
        grid on
        hold off
        
    end
    
end


