% TestLinfitErrorBoundsOnYhat.m   Feb. 2006 L Braile, revised 3/3/2011
% Test and plot the linfit least squares straight
% line function with sample data; also calculate and plot the
% 99% confidence or error bounds on the predicted values of
% y (Y-hat); these bounds can be used to objectively
% identify outliers.
ArcAgeSlabDipData    % use arc Age data
x = A(:,1);  % First arc age (duration of subduction) in millions of years
y = A(:,2);  % Second column is slab dip angle in degrees
mode = 0;
npts = length(x);
sigmay = ones(1,npts); % assume uncertainties in y are all = 1
% calculate and print least squares coefficients, sigmas and r
[a,sigmaa,b,sigmab,r] = linfit(x,y,sigmay,npts,mode) % use linfit
% to calculate the least squares coefficients and correlation
% coefficient
% Calculate the confidence or error bounds on Y-hat
yt = a + b*x;  % calculate the theoretical values of y at the data points
SSE = sum((y-yt).^2);
MSE = SSE/(npts-2); % note d.f. is npts - 2; 
s = sqrt(MSE);  % the MSE is used as an estimate of sigma^2
tcoeff = 3; % ~ corresponds to 99% confidence (should be 2.6 for 
% infinite size sample, ~ 2.75 for npts = 30; often rounded to 3)
%tcoeff = 2; % ~ corresponds to 95% confidence (should be 1.96 for 
% infinite size sample, ~ 2.04 for npts = 30; often rounded to 2)
%tcoeff = 2.2;  % value from tables for 95% and npts = 13 (core data)
xt = [0:1:250]; % calculate x locations for plotting
yhat = a + b*xt; % calculate y-hat at x locations for plotting
% CSSX = sum(((x - mean(x)).^2)); % These two equations for CSSX
CSSX = sum(x.^2) - ((sum(x)).^2)/npts; % are equivalent
bound = tcoeff * s * sqrt(((1/npts) + ((xt - mean(x)).^2)/CSSX)); % calculate
% bounds on y-hat
% plot the data and least squares line
plot(x,y,'ro','markersize',12,'linewidth',3) % Plot the data
set(gca,'fontsize',16,'linewidth',2)
hold on
% plot bounds
plot(xt, yhat + bound,'--b')
plot(xt, yhat - bound,'--b')
xlabel('Arc Age (Duration of Subduction, Ma)','fontsize',16)
ylabel('Slab Dip Angle (Degrees)','fontsize',16)
title('Arc Age vs. Slab Dip Angle (Data from Jarrard, 1986)','fontsize',16)
yt = feval('sline',x,a,b); % Calculate the fitted line
plot(x,yt,'k-','linewidth',3) % Plot the fitted line
axis([0 250 0 80]); % Set axis limits for plot Arc age data
hold off