;% lsqfitgm.m                                     by:  Edward T Peltzer, MBARI
;%                                                revised:  2000 Jan 31.
;%
;% 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)
 
;function [m,b,r,sm,sb]=lsqfitgm(X,Y)
PRO LSQFITGM,x,y,m,b,r,sm,sb,yfit

;% Determine slope of Y-on-X regression

;[my] = lsqfity(X,Y)
res1=LINFIT(x,y,sigma=sigma)
my=res1[1]

;% Determine slope of X-on-Y regression

;[mxi] = lsqfitxi(X,Y)
res2=LINFIT(y,x,sigma=sigma2)
mxi=1./res2[1]

;% Calculate geometric mean slope

m = SQRT(my * mxi)

IF (my LT 0) AND (mxi LT 0) THEN m=-m

;% Determine the size of the vector
n = N_ELEMENTS(x)

;% Calculate TOTALs and means
Sx = TOTAL(X)
Sy = TOTAL(Y)
xbar = MEAN(x)
ybar = MEAN(y)

;% Calculate geometric mean intercept
b=ybar-m*xbar

;% Calculate more TOTALs
Sxy = (x*y)
Sx2 = TOTAL(x^2)
Sy2 = TOTAL(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/mxi)

IF (my LT 0 AND mxi LT 0) THEN r = -r

diff= Y - b - m * X

s2 = TOTAL(diff * diff)/(n-2)
sm = SQRT(n*s2/den)
sb = SQRT(Sx2*s2/den)
yfit=m*x+b

;stop
END