% lsqfitma.m                                     by:  Edward T Peltzer, MBARI
%                                                revised:  2007 Sep 05.
% 
% M-file to calculate a "MODEL-2" least squares fit.
%
%     The line is fit by MINIMIZING the NORMAL deviates.
%
%     The equation of the line is:     y = mx + b.
%
%     This line is called the MAJOR AXIS.  All points are given EQUAL
%       weight.  The units and range for X and Y must be the same.
%     Equations are from York (1966) Canad. J. Phys. 44: 1079-1086;
%       re-written from Kermack & Haldane (1950) Biometrika 37: 30-41;
%       after a derivation by Pearson (1901) Phil. Mag. V2(6): 559-572.
%
%     Data are input and output as follows:
%
%	    [m,b,r,sm,sb] = lsqfitma(X,Y)
%
%             X    =    x data (vector)
%             Y    =    y data (vector)
%
%             m    =    slope
%             b    =    y-intercept
%             r    =    correlation coefficient
%             sm   =    standard deviation of the slope
%             sb   =    standard deviation of the y-intercept
%
%     Note that the equation passes through the centroid:  (x-mean, y-mean)
 
function [m,b,r,sm,sb]=lsqfitma(X,Y)

% Determine the size of the vector
 
n = length(X);
 
% Calculate sums and other re-used expressions
 
Sx = sum(X);
Sy = sum(Y);
xbar = Sx/n;
ybar = Sy/n;
U = X - xbar;
V = Y - ybar;
 
Suv = sum(U .* V);
Su2 = sum(U .^2);
Sv2 = sum(V .^2);
 
sigx = sqrt(Su2/(n-1));
sigy = sqrt(Sv2/(n-1));
 
% Calculate m, b, r, sm, and sb
 
m = (Sv2 - Su2 + sqrt(((Sv2 - Su2)^2) + (4 * Suv^2)))/(2 * Suv);
b = ybar - m * xbar;
r = Suv / sqrt(Su2 * Sv2);
 
sm = (m/r) * sqrt((1 - r^2)/n);
sb1 = (sigy - sigx * m)^2;
sb2 = (2 * sigx * sigy) + ((xbar^2 * m * (1 + r))/r^2);
sb = sqrt((sb1 + ((1 - r) * m * sb2))/n);