%% York Regression
% A bivariate weighted regression, useful when x and y errors are known
% from York 2004
% see https://amt.copernicus.org/articles/11/1233/2018/amt-11-1233-2018.pdf
% can be calculated using origin as well


function [y_a err_y_a y_b err_y_b]=f_york(x,y, xerr, yerr) % given the inputs 
% function [y_b err_y_b]=f_york(x,y, xerr, yerr) % given the inputs 

xbar = mean(x); % mean of x
ybar = mean(y); % mean of y
n = length(x); % number of data points
tcrit = 1.99; % t critical value, changes with n

% Regression calculations

% calculate x and y omegas
omega_x = (1./(xerr/2).^2);
omega_y = (1./(yerr/2).^2);

% r is the correlation coeeficient between x and y errors
y_r_m = corrcoef(xerr,yerr);
y_r = y_r_m(1,2);
% b is the univariate regression slope (ordinary least squares)
ols = fitlm(x,y,'linear');
ols_b = table2array(ols.Coefficients(2,1)); % slope is in coeffcients table
% breakdown into smaller components
y_W = (omega_x.*omega_y)./(omega_x + ols_b.^2.*omega_y - 2*ols_b.*y_r*sqrt(omega_x.*omega_y));
wxbar = sum(y_W.*x)/sum(y_W);
wybar = sum(y_W.*y)/sum(y_W);
y_U = x - wxbar;
y_V = y - wybar;
y_Beta = (y_W.*((y_U./omega_y)+(ols_b.*y_V./omega_x)-((ols_b.*y_U+y_V).*(y_r./sqrt(omega_x.*omega_y)))));

% calculate slope and intercept
y_b = (sum(y_W.*y_Beta.*y_V))/(sum(y_W.*y_Beta.*y_U)); % slope
y_a = ybar - y_b*xbar; % intercept

% Calculate uncertainties
y_x = xbar + y_Beta;
y_xbar = mean(y_x);
y_u = y_x - y_xbar;
% 2SD using standard errors
err_y_b = tcrit*sqrt((1./sum(y_W.*y_u.^2)));
err_y_a = tcrit*sqrt(((1./sum(y_W))+y_xbar^2*(err_y_b/tcrit)^2));

% y_T = table(y_a, err_y_a, y_b, err_y_b);
% y_T.Properties.VariableNames = {'intercept','intercept 2SD','slope','slope 2SD'}
end
