function [a, yest, sst, rr, F] = xyplot(x,y); % from Middleton, p. 60
% with modifications (mostly comments) by L Braile, Feb. 2005
% [a, yest, sst, rr] = xyplot(x,y) revised Feb. 2006
% plots y vs x and calculates linear regression, correlation 
% and ANOVA
%   x vector of independent variable
%   y vector of dependent variable
%   a - vector of fit coefficients: 
%       a(1) is slope, a(2) is y-intercept
%   yest - vector of estimated fit points
%   sst  - total sums of square
%   rr -  square of correlation coefficient
%   F  - calculate F statistic
N = length(x);
Z = [x y];
r = corrcoef(Z);		% calculate rr, r(1,2) is the corr. coeff. 
rr = r(1,2)*r(1,2);     % rr is r^2, the "goodness of fit"
sdy = std(y);			% calculate total sums of squares
sst = sdy*sdy*(N-1);
clc; % clears the display window
% The following code illustrates how to make an output data table and
% how to include the value of calculated variables in the output.
fprintf('Analysis of Variance Table.\n\n');
fprintf('Source       d.f.     Sums of Squares  Mean Squares\n');
fprintf('---------------------------------------------------\n');
fprintf('Regression    1       %8.2f         %8.2f\n', (rr*sst),(rr*sst) );
fprintf('Residuals    %2.0f       %8.2f         %8.2f\n',N-2,(1-rr)*sst, (1-rr)*sst/(N-2) );
fprintf('---------------------------------------------------\n');
fprintf('Total        %2.0f       %8.2f\n\n',N-1, sst);
m = 1;
a = polyfit(x,y,m);					 % calculate regression
yest = polyval(a,x);
F = (N-2)*rr/(1-rr);
fprintf('F = %g\n', F);
fprintf('Correlation coeff: r = %g\n',r(1,2) );
fprintf('Equation: y = %g + %g x\n\n',a(2),a(1) );
%fprintf('Click on Window to see figure\n');
plot(x,y,'o',x,yest,'-');
xlabel('x'), ylabel('y'), title('Linear Regression');