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 25 0 20]); 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 plot(x,yhat,'-r') 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 25 0 20]); 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 25 0 20]); 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 25 0 20]); 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