function [ ] = linreg( )
% RESET RANDOM NUMBER GENERATOR
rng('default');
% GENERATE DATA WITH LINEAR TREND
r1=10;
r2=10;
ls=0;
le=10;
n=101;
x=linspace(ls,le,n);
x=x+r1*rand(size(x));
y=linspace(ls,le,n);
y=y+r2*rand(size(y));

%CLEAR ALL GRAPHS
clf;
%PLOT 1 Y-on-X FIT
subplot(2,2,1)  % 2 x 2 arrangement of plots: 1 - left upper
title('Data and Y-on-X Fit');
hold on
axis( [0 20 0 20],'square');
p=polyfit(x,y,1);   %Y on X regression, 1 for linear
yfit=p(1)*x +p(2);
yresid=y-yfit;
plot(x,y,'ko',x,yfit,'-g')
pause
for i=1:n
    xp(1)=x(i);
    xp(2)=x(i);
    yp(1)=y(i);
    yp(2)=yfit(i);
    plot(xp,yp,'-g');
    hold on
end
pause
% test alternative computation of Y on X regression
meanX=mean(x);
meanY=mean(y);
SXX = sum((x-meanX).^2);
SYY = sum((y-meanY).^2);
%SXY = sum( (x-mean(x)).*(y-mean(y) ) );
SXY = sum( (x-meanX).*(y-meanY ) );
slope = SXY/SXX
intercept = meanY - slope*meanX
yhat = slope*x + intercept
 

hold on
SSresidYonX=sum(yresid.^2)
SStotalY=(length(y)-1)*var(y)
rsqYonX=1-SSresidYonX/SStotalY
pause

%PLOT 2 X-on-Y FIT
subplot(2,2,2)
title('Data and X-on-Y Fit')
hold on
axis( [0 20 0 20],'square');

hold on
q=polyfit(y,x,1);
xfit=q(1)*y+q(2);
xresid=x-xfit;
plot(x,y,'ko',xfit,y,'-r')
pause
for i=1:n
    xp(1)=x(i);
    xp(2)=xfit(i);
    yp(1)=y(i);
    yp(2)=y(i);
    plot(xp,yp,'-r');
    hold on
end
pause
%plot(x,y,'ko',xfit,y,'-r')
hold on
SSresidXonY=sum(xresid.^2)
SStotalX=(length(x)-1)*var(x)
rsqXonY=1-SSresidXonY/SStotalX
pause

%PLOT 3 BISECTOR OF Y-on-X and X-on-Y FITS
subplot(2,2,3)
hold on
title('Data and Bisector Fit');
hold on
axis([0 20 0 20],'square');
hold on
xbar=mean(x);
ybar=mean(y);
theta1=atan(p(1));
theta2=atan(1/q(1));
theta=(theta1+theta2)/2;
b=tan(theta);
xs=linspace(ls,le+r1,n);
ys=b*(xs-xbar)+ybar;
SSresidBis=0.;
plot(x,y,'ko',xs,ys,'-b')
pause

for i=1:n
    % solve (y -y(i )=-1/b*(x-x(i)) and
    %       (ys-ybar)=b   *(x-xbar) for x and y
    % to get intersection
    xp(1)=x(i);
    yp(1)=y(i);
    xp(2)=(b*b*xbar-b*ybar+y(i)*b+x(i))/(1+b*b);
    yp(2)=(-b*xbar+ybar+y(i)*b*b+b*x(i))/(1+b*b);
    SSresidBis=SSresidBis+(xp(1)-xp(2))^2+(yp(1)-yp(2))^2;
    plot(xp,yp,'-b')
    hold on
end
pause

SSresidBis
pause

%PLOT 4 ALL FITS
subplot(2,2,4)
title('Data and Fits')
hold on
axis([0 20 0 20],'square');
hold on
plot(x,y,'ko',x,yfit,'-g',xfit,y,'-r',xs,ys,'-b');
hold on
hl=legend('Data','Y-on-X','X-on-Y','Bisector');
set(hl,'Location','SouthEast');
hold on;

end