%
% 2D Linear Regression Fit of y = a * x + b
%
%      [ p, h ] = linreg( y, x, w )
%
% where  p = fitted parameters = [ b a ]
%        h = handles of the five plots:
%            h(1) = data and regression line
%            h(2) = distribution of residuals
%            h(3) = residuals against x
%            h(4) = residuals against y
%            h(5) = parameter confidence region
%        y = dependent data vector
%        x = independent data vector
%        w = data weights vector
%
% the size of y, x, and w should be the same
%
% Bob Newell, February 1996
%
%------------------------------------------------------
%
function [ p, h ] = linreg( y, x, w )
%
fontsz = 16;
%
disp( 'LINEAR REGRESSION' );
h = zeros( 5, 1 );
% regression matrix
x = x(:);
X = [ ones( size(x) ) x ];
% data vector
y = y(:);
% some weighting to try and compensate
% for the transformation
if length( w ) == length( y ),
  W = diag( w );
else,
  W = diag( ones( size( y ) ) );
end;
% perform the linear regression
p = inv( X' * W * X ) * X' * W * y;
%---------------------------------------------------
% the answers
%
disp( 'Fitted Parameters:' );
disp( [ '    Intercept = ' num2str( p(1) ) ] );
disp( [ '    Slope     = ' num2str( p(2) ) ] );
h(1) = figure;
plot( x, y, '+' )
set( gca, 'fontsize', fontsz );
hold on
ax = axis;
plot( [ ax(1:2) ], [ p(2) * ax(1:2) + p(1) ] )
xlabel( 'X Data' )
ylabel( 'Y Data' )
hold off
%----------------------------------------------------
% residual plots
%
resid = y - X * p;
% first the distribution plot
h(2) = resdist( resid );
set( gca, 'fontsize', fontsz );
xlabel( 'Residuals' )
ylabel( 'Frequency' )
% plot against x data
h(3) = figure;
plot( X(:,2), resid, '*' )
set( gca, 'fontsize', fontsz );
ax = axis;
line( [ ax(1:2) ], [ 0 0 ] )
xlabel( 'X Data' )
ylabel( 'Residuals' )
% plot against y data
h(4) = figure;
plot( y, resid, '*' )
set( gca, 'fontsize', fontsz );
ax = axis;
line( [ ax(1:2) ], [ 0 0 ] )
xlabel( 'Y Data' )
ylabel( 'Residuals' )
%-------------------------------------------------
% calculate percentage explained variation
%
% total variation
sostot = ( y - mean(y) )' * ( y - mean(y) );
% residual variation
sosleft = resid' * resid;
% percentage explained
pcfit = 100 * ( sostot - sosleft ) / sostot;
disp( [ 'Percent Variance Explained by Fit = ' ...
             num2str( floor( pcfit ) ) ] );
%--------------------------------------------------
% correlation of parameters
%
% unexplained variance (estimate of error variance,
% but of course includes model mismatch)
var = resid' * resid / ( length( resid ) - 2 );
% parameter covariance matrix
parvar = inv( X' * X ) * var;
% parameter correlation coefficient (normalised)
d = diag( parvar );
corr = parvar ./ sqrt( d * d' );
disp( [ 'Correlation between Parameters = ' ...
                           num2str( corr(1,2) ) ] );
%-----------------------------------------------------
% parameter confidence limits
%
% read from tables, for 95% and dof=96, t=1.99
%
cl = 1.99 * sqrt( d );
disp( 'Confidence Limits:' )
disp( [ '    Intercept = ' num2str( cl(1) ) ] )
disp( [ '    Slope     = ' num2str( cl(2) ) ] )
%-----------------------------------------------------
% parameter confidence region
%
% read from tables, for 95% and dof=(2,96), F=7.53
%
% sos on the 95% contour
sos95 = sosleft * ( 1 + 2 * 7.53 / 96 );
% plot these contour values
crv = [ 1:4 ] * sos95;
% evaluate sos at these parameter grid points
crx = [ 1:21 ] * cl(1) / 2 + p(1) - 5 * cl(1);
cry = [ 1:21 ] * cl(2) / 2 + p(2) - 5 * cl(2);
% calculate the sos values
s = zeros( 21, 21 );
for i = 1:21,
  for j = 1:21,
    r = y - X * [ crx(i) cry(j) ]';
	s(i,j) = r' * r;
  end;
end;
% then do the plot of contours
h(5) = figure;
contour( crx, cry, s, crv )
set( gca, 'fontsize', fontsz );
grid
xlabel( 'Intercept' )
ylabel( 'Slope' )
% add blob for the parameter estimates
hold on
ax = axis;
blobx = 0.020 * ( ax(2) - ax(1) );
bloby = 0.025 * ( ax(4) - ax(3) );
xblob = blobx * [ 1 -1 -1 1 ]' + p(1);
yblob = bloby * [ -1 -1 1 1 ]' + p(2);
fill( xblob, yblob, 'k' )
% add square for confidence limits
line( [p(1)-cl(1) p(1)-cl(1)], [p(2)-cl(2) p(2)+cl(2)] );
line( [p(1)+cl(1) p(1)+cl(1)], [p(2)-cl(2) p(2)+cl(2)] );
line( [p(1)+cl(1) p(1)-cl(1)], [p(2)+cl(2) p(2)+cl(2)] );
line( [p(1)+cl(1) p(1)-cl(1)], [p(2)-cl(2) p(2)-cl(2)] );
hold off
return;
%-------------------------------------------------
% the end
