% lsqbisec.m by: Edward T Peltzer, MBARI
% revised: 2016 Mar 17.
%
% M-file to calculate a "MODEL-2" least squares fit.
%
% The SLOPE of the line is determined by calculating the slope of the line
% that bisects the minor angle between the regression of Y-on-X and X-on-Y.
%
% The equation of the line is: y = mx + b.
%
% This line is called the LEAST SQUARES BISECTOR.
%
% See: Sprent and Dolby (1980). The Geometric Mean Functional Relationship.
% Biometrics 36: 547-550, for the rationale behind this regression.
%
% Sprent and Dolby (1980) did not present a statistical treatment for the
% estimation of the uncertainty limits for the least squares bisector
% slope, or intercept.
%
% 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] = lsqbisec(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]=lsqbisec(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 the least squares bisector slope
theta = (atan(my) + atan(mx)) / 2;
m = tan(theta);
% 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 the least squares bisector intercept
b = ybar - m * xbar;
% Calculate more sums
Sxy = sum(X .* Y);
Sx2 = sum(X.^2);
% Calculate re-used expressions
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);