% lsqfitgm.m                                     by:  Edward T Peltzer, MBARI
%                                                revised:  2007 Apr 28.
% 
% M-file to calculate a "MODEL-2" least squares fit.
%
%     The SLOPE of the line is determined by calculating the GEOMETRIC MEAN
%       of the slopes from the regression of Y-on-X and X-on-Y.
%
%     The equation of the line is:     y = mx + b.
%
%     This line is called the GEOMETRIC MEAN or the REDUCED MAJOR AXIS.
%
%     See Ricker (1973) Linear regressions in Fishery Research, J. Fish.
%       Res. Board Can. 30: 409-434, for the derivation of the geometric
%       mean regression.
%
%     Since no statistical treatment exists for the estimation of the
%       asymmetrical uncertainty limits for the geometric mean slope,
%       I have used the symmetrical limits for a model I regression
%       following Ricker's (1973) treatment.  For ease of computation,
%       equations from Bevington and Robinson (1992) "Data Reduction and
%       Error Analysis for the Physical Sciences, 2nd Ed."  pp: 104, and
%       108-109, were used to calculate the symmetrical limits: sm and sb.
%
%     Data are input and output as follows:
%
%	    [m,b,r,sm,sb] = lsqfitgm(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)
%
%     WARNING:  Both lsqfitx.m and lsqfity.m must be present in a directory
%               that is in your MATLABPATH in order for this algorithm to
%               execute properly.

function [m,b,r,sm,sb]=lsqfitgm(X,Y)

% Determine slope of Y-on-X regression

[my] = lsqfity(X,Y);

% Determine slope of X-on-Y regression

[mx] = lsqfitx(X,Y);

% Calculate geometric mean slope

m = sqrt(my * mx);

if (my < 0) & (mx < 0)
	m = -m;
end

% Determine the size of the vector
 
n = length(X);
 
% Calculate sums and means
 
Sx = sum(X);
Sy = sum(Y);
xbar = Sx/n;
ybar = Sy/n;

% Calculate geometric mean intercept

b = ybar - m * xbar;

% Calculate more sums

Sxy = sum(X .* Y);
Sx2 = sum(X.^2);
Sy2 = sum(Y.^2);

% Calculate re-used expressions

num = n * Sxy - Sx * Sy;
den = n * Sx2 - Sx^2;

% Calculate r, sm, sb and s2

r = sqrt(my / mx);

if (my < 0) & (mx < 0)
	r = -r;
end

diff = Y - b - m .* X;

s2 = sum(diff .* diff) / (n-2);
sm = sqrt(n * s2 / den);
sb = sqrt(Sx2 * s2 / den);