function [S,EVR,MSE]=regresion(X,Y,method)
% 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
% from http://www.mathworks.com/matlabcentral/newsreader/view_thread/158607
% Originally created under /Users/Mau/Documents/FALL2010/DataAnalysis_Gille/Project3/QscatvsWindsat


% 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]=corr(x,y);

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

yhat = x*S; % model
EVR = var(yhat)/var(y); % "Explained" Variance Ratio
plot(x,yhat,'r');
% e = y-yhat; % error
% MSE = mse(e); % mean-squared error

plot(x,y,'ok','markersize',6,'markerfacecolor','k');
set(gca,'fontsize',16);
hold on
plot(x,yhat,'r','linewidth',2);
legend(['Correlation: r = ', num2str(r, '%6.2f'),', pval = ', num2str(pval, '%6.4f')],...
    [' y =',num2str(S, '%6.2f'),' x,  MSE = ',num2str(MSE,'%6.2f' )],...
    'location','NE','fontsize',16);
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);
% figure, plot(x,res,'+')
% title('Residuals for the Linear Model')


plot(x,y,'ok','markersize',6,'markerfacecolor','k');
set(gca,'fontsize',16);
hold on
plot(x,y_fit,'r','linewidth',2);
legend(['Correlation: r = ', num2str(r, '%6.2f'),', pval = ', num2str(pval, '%6.4f')],...
    [' y =',num2str(p(1), '%6.2f'),' x + ',num2str(p(2), '%6.2f'), ', MSE = ',num2str(MSE,'%6.2f' )],...
    'location','NE','fontsize',16);
grid on
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==2
% No intercept linear model (least squares solution)
[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,y,'o','markersize',10,'markerfacecolor',[0 0.5 0.6],'color',[0 0.5 0.6]);
set(gca,'fontsize',20);
hold on
plot(x,yhat,'k','linewidth',2);
hold on
grid on
    
% Linear regression (least squares solution)   
% Calculate fit parameters
[p,ErrorEst] = polyfit(x,y,1);
% 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);


plot(x,y_fit,'color',[0.5 0.5 0.5],'linewidth',2);
legend(['Correlation: r = ', num2str(r, '%6.2f'),', pval = ', num2str(pval, '%6.4f')],...
    [' y =',num2str(S, '%6.2f'),' x, MSE = ',num2str(MSE,'%6.2f' )],...
    [' y =',num2str(p(1), '%6.2f'),' x + ',num2str(p(2), '%6.2f'), ' , MSE = ',num2str(MSE1,'%6.2f' )],...
    'location','NE','fontsize',16);
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)

end


