function linear_regression_YORK
%
% linear_regression_YORK is a function that fits a straight line 
%   y = a + b*x  to a set of X,Y data points using York's Linear Iterative
%   equation.  It is assumed that the measurements have estimates of standard 
%   deviations sX and sY and covariance sXY.
%   The function reads an ASCII textfile containing coordinate pairs (X,Y)
%   standard deviations (sX,sY) and covariance (sXY).  
%   In addition point number is included in the textfile.
%   Part of an ASCII textfile is shown below
%
% Neri_Line_Data_2.txt
% Line data from Neri et al (1989) but with covariances
%
% % serial  X     sX            Y     sY            sXY
%     1    0.0   0.031622777   5.9   1.0           0.031274926
%     2    0.9   0.031622777   5.4   0.745355993  -0.020506097
%     3    1.8   0.044721360   4.4   0.5          -0.004986432
%     4    2.6   0.035355339   4.6   0.353553391   0.001237500
%     5    3.3   0.070710678   3.5   0.223606798  -0.000901249
%     6    4.4   0.111803399   3.7   0.223606798  -0.016500000
%     7    5.2   0.129099445   2.8   0.119522861   0.000339467
%     8    6.1   0.223606798   2.8   0.119522861   0.019804058
%     9    6.5   0.745355003   2.4   0.1          -0.024969426
%    10    7.4   1.0           1.5   0.044721360  -0.000044721
% % end of data
%
%   Results are written to a textfile having the same path and name as the 
%   data file but with the extension ".out"

%===============================================================================
% Function:  linear_regression_YORK
%
% Author:
%  Rod Deakin, 
%  BONBEACH VIC 3196
%  AUSTRALIA
%  email: randm.deakin@gmail.com
%
% Date:
%  Version 1.0  27 June 2017
%  Version 1.1  25 July 2017
%
% Remarks:
%  This function reads numeric data from a textfile containing coordinate
%  pairs (Y,Y) standard deviation and covariances associated with each pair 
%  and computes the parameters a and b of a line of best fit y = a + b*x 
%  using York's Linear Iterative equation. 
%  Results are written to a textfile having the same path and name as the 
%  data file but with the extension ".out"
%
% Arrays:
%  alpha     - (n,1) vector of square-roots of products of weights wX, wY
%  beta      - (n,1) vector
%  rXY       - (n,1) vector of correlations
%  serial    - (n,1) vector of station serial numbers
%  sX        - (n,1) vector of st. dev's of X coordinates
%  sXY       - (n,1) vector of covariances between Y and Y coordinates
%  sY        - (n,1) vector of st. dev's of Y coordinates
%  U         - (n,1) vector of U centroidal coordinates (U = X - Xc)
%  V         - (n,1) vector of V centroidal coordinates (V = Y - Yc)
%  vX        - (n,1) vector of X residuals
%  vY        - (n,1) vector of Y residuals
%  wX        - (n,1) vector of weights wX = 1/sX^2
%  wY        - (n,1) vector of weights wY = 1/sY^2
%  X         - (n,1) vector of measured X coordinates 
%  Y         - (n,1) vector of measured Y coordinates
%  x         - (1,n) vector of computed x coordinates
%  y         - (1,n) vector of computed y coordinates
%  xx,yy     - vectors of coordinates for plotting Line of Best Fit
%  W         - (n,1) vector of weight coefficients W(k)
%
% Variables:
%  a         - point on y-axis where straight line intersects
%  b         - b = tan(theta) = (Y2-Y1)/(X2-X1) = slope of straight line
%  b0        - approximate value of slope of straight line
%  beta_bar  - factor in computation of variance estimates of b and a
%  D         - denominator in calculation of variance of gradient
%  eps       - small quantity for testing for convergence
%  iter      - iteration count
%  maxX      - maximum X coordinate
%  minX      - minimum X coordinate
%  maxY      - maximum Y coordinate
%  minY      - minimum Y coordinate
%  n         - number of measurement pairs
%  k         - integer counter
%  S         - sum of weighted squared residuals
%  sa        - st.dev of gradient a
%  sb        - st.dev of gradient b
%  s2b       - estimate of variance of b
%  s2a       - estimate of variance of a
%  Xc,Yc     - X and Y centroids
%
% References:
%  York, D., 1968, 'Least-squares fitting of a straight line with 
%      correlated errors', Earth and Planetary Science Letters (1969),
%      pp. 320-324, North-Holland Publishing Company, Amsterdam.
%
%===============================================================================

%-------------------------------------------------------------------------------
% 1. Call the User Interface (UI) to choose the input data file name
% 2. Concatenate strings to give the path and file name of the input file
% 3. Strip off the extension from the file name to give the rootName
% 4. Add extension ".out" to rootName to give the output filename
% 5. Concatenate strings to give the path and file name of the output file
%-------------------------------------------------------------------------------
filepath = strcat('d:\temp\','*.txt');
[infilename,inpathname] = uigetfile(filepath);
infilepath  = strcat(inpathname,infilename);
rootName    = strtok(infilename,'.');
outfilename = strcat(rootName,'.out');
outfilepath = strcat(inpathname,outfilename);

%-------------------------------------------------------------------------------
% 1. Load the data into an array whose name is the rootName
% 2. Copy columns of data into individual arrays
%-------------------------------------------------------------------------------
fileTemp = load(infilepath);
serial  = fileTemp(:,1);
X   = fileTemp(:,2);
sX  = fileTemp(:,3);
Y   = fileTemp(:,4);
sY  = fileTemp(:,5);
sXY = fileTemp(:,6);
%-------------------------------------------------------------------------------
% Determine the number of measurement pairs n and form vectors of weights,
% correlations and square roots of products of weights
%-------------------------------------------------------------------------------
n     = length(X);
wX    = zeros(n,1);
wY    = zeros(n,1);
rXY   = zeros(n,1);
wX    = 1./sX.^2;
wY    = 1./sY.^2;
rXY   = sXY./(sX.*sY);
alpha = sqrt(wX.*wY);

%-------------------------------------------------------------------------------
% Print heading and given data
%-------------------------------------------------------------------------------
% open the output file print the data
fidout  = fopen(outfilepath,'wt');
fprintf(fidout,'\n\n Linear Regression Using YORK''S Linear Iterative equation');
fprintf(fidout,'\n ========================================================');
fprintf(fidout,'\n solution for the parameters a and b in the equation for the line of best fit');
fprintf(fidout,'\n y = a + b*x  where b = tan(theta) is the gradient and');
fprintf(fidout,'\n a is the intercept on the y-axis');
fprintf(fidout,'\n\n Input Data');
fprintf(fidout,'\n Serial      X       sX              Y       sY            sXY');
for k = 1:n
    fprintf(fidout,'\n %4d %12.4f  %8.6f  %12.4f  %8.6f     %9.6f',...
            serial(k,1),X(k,1),sX(k,1),Y(k,1),...
            sY(k,1),sXY(k,1));
end
fprintf(fidout,'\n\n there are n = %2d measurement pairs (the X,Y coordinates),',n); 
fprintf(fidout,'\n sX,sY are estimates of measurement standard deviations and'); 
fprintf(fidout,'\n sXY are covariances'); 

%-------------------------------------------------------------------------------
% Calculate approximate solution for slope b = tan(theta) and intercept a.  
% assuming all observations are of unit weight and uncorrelated
%-------------------------------------------------------------------------------
% calculate the X,Y coords of the centroid
Xc = sum(X)/n;
Yc = sum(Y)/n;
% calculate the centroidal coordinates U,V
U = X-Xc;
V = Y-Yc;
% calculate sums of squares and products of U,V coordinates
U2 = sum(U.^2);
V2 = sum(V.^2);
UV = sum(U.*V);
% calculate approximation to b
b0 = ((V2-U2) + sqrt((V2-U2)^2 + 4*(UV)^2)) / (2*UV);

%-------------------------------------------------------------------------------
% calculate gadient b using iterative linear equation of York
%-------------------------------------------------------------------------------
b    = b0;
iter = 1;
eps  = 1e-9;
while 1
  W = (alpha.^2)./(b^2.*wY + wX - 2*b.*rXY.*alpha);
  Xc = sum(W.*X)/sum(W);  
  Yc = sum(W.*Y)/sum(W);
  U = X - Xc;
  V = Y - Yc; 
  beta = W.*(U./wY + b.*V./wX - (rXY./alpha).*(b.*U + V));
  b_next = sum(W.*beta.*V)/sum(W.*beta.*U);
  if abs(b_next - b) < eps
    break;
  end;
  iter = iter+1;
  b = b_next;
  if iter > 20
    fprintf('Iteration for gradient b failed to converge after 20 iterations');
    break;
  end;
end;
% set final values of b and a
b = b_next;
a = Yc - b*Xc;
% compute adjusted coords
x = beta + Xc;
y = b.*beta + Yc;
% compute residuals
vX = x - X;
vY = y - Y;
% compute sum of squares
S = sum(W.*(b.*U - V).^2);
% compute beta_bar and D
beta_bar = sum(W.*beta)/sum(W);
D = 1/b*sum(W.*U.*V-W.^2.*rXY./alpha.*(b.*U-V).^2) + 4*sum(W.*(beta-U).*(beta-beta_bar));
% compute variance of b
s2b = 1/D^2*sum(W.^2.*(U.^2./wY - 2.*rXY./alpha.*U.*V + V.^2./wX));
% compute variance of a
s2a = (2*beta_bar+Xc)^2*s2b + 2/D*(2*beta_bar+Xc)*beta_bar + 1/sum(W);

%-------------------------------------------------------------------------------
% Print solutions
%-------------------------------------------------------------------------------
fprintf(fidout,'\n\n Solutions after %2d iterations',iter);
fprintf(fidout,'\n =============================');
fprintf(fidout,'\n slope of line b = tan(theta) = % 12.9f',b');
fprintf(fidout,'\n y-axis intercept           a = % 12.9f',a');

fprintf(fidout,'\n\n Vector of residuals');
fprintf(fidout,'\n   k    vX             vY');
for k = 1:n
    fprintf(fidout,'\n %3d  % 12.9f   % 12.9f',k, vX(k,1),vY(k,1));
end

fprintf(fidout,'\n\n sum of weighted squared residuals S = %15.9f',S);
fprintf(fidout,'\n\n variance of gradient b    = %15.9f',s2b);
fprintf(fidout,'\n variance of y-intercept a = %15.9f',s2a);
fprintf(fidout,'\n\n Reference: York, D., 1968, ''Least-squares fitting of a straight line with'); 
fprintf(fidout,'\n                correlated errors'', Earth and Planetary Science Letters (1969),');
fprintf(fidout,'\n                pp. 320-324, North-Holland Publishing Company, Amsterdam.');

fprintf(fidout,'\n\n');
 
% close the output file
fclose(fidout);

%-------------------------------------------------------------------------------
% plot data points and line of best fit
%-------------------------------------------------------------------------------
% set max, min of measured coordinates
maxX   = max(X);
minX   = min(X);
maxY   = max(Y);
minY   = min(Y);
xx = linspace(minX,maxX,500);
yy = (b.*xx) + a;

% Select Figure window and clear figure
figure(1);
clf(1);
hold on;
grid on;
box on;
axis ([minX,maxX,minY,maxY],'equal');

% plot line of best fit and then the data points with a star (*)
plot(xx,yy,'-k');
plot(X,Y,'*k');
% plot the centroid
plot(Xc,Yc,'ok');

% anotate the plot
title('Line of Best Fit: measurement(*), centroid(o)'); 
xlabel('X measurements');
ylabel('Y measurements');
