%%% jcaleyIII.m %%%
%%% written by Patrick F. Dunn %%%
%%% copyright 2003 %%%
%%% modified on 2/15/2005 and later %%%

% this m-file reads a three-column calibration data file of x, y, y-error,
% performs a linear regression analysis to obtain the calibration fit,
% then plots various confidence intervals associated with the fit and
% with the use of the fit. It also detemines the percentage uncertainty in the
% x-from-y estimate for the range of y values and plots the percent uncertainty
% in the x-from-y estimate versus x. This percentage uncertainty becomes
% infinite when x=0, so the percentage uncertainty at x=0 is not plotted.
% Reference: Chpt.12: John Mandel's "The Statistical Analysis of Experimental 
% Data"

% if casenox=1 then confidence intervals are determined and plotted for:
%
% case 1: uncertainty in y due to its measurement error + fit error;
%         for situation when fit is NOT used to to find x values
%         (confidence intervals are two lines)
%         This situation also models the case of controlled x (Berkson case)
%         which assumes both x and y errors that are statistically independent
%         although interpretation of the confidence limits is different 
%         (see p.295)
% case 2: uncertainty in x determined by fit and average y value of n 
%         replications;
%         for situation when fit is used to find x value
%         that is determined from n replications of y
%         (confidence intervals are two hyperbolas)
%			 a ypick value is user specified; uncertainty in 
%         xfromypick is estimated and shown on plot

clear all; clf;

caseno1=1; %if =1 this will plot the confidence intervals for case 1
caseno2=1; %if =1 this will plot the confidence intervals for case 2
% cases 1 and/or 2 can be plotted

% JLO
uncx = [];
uncy = [];

% the input data file is read
filename = input('enter the filename w/o its extension: ','s');
ext = input('enter the files extension, e.g., dat: ','s');
% filename = 'data';
% ext = 'txt';
eval(['load ',filename,'.',ext]);
x = eval([filename,'(:,1)']);
y = eval([filename,'(:,2)']);
ey = eval([filename,'(:,3)']);
xes=size(x);
% JLO y2 not used to scale uncertainty axis x2=x(2);

% the percent confidence is set
% typically 95%
P = input('enter the percent confidence: P(%) = ');

% the x,y data is fit and evaluated using least-squares linear regression
p = polyfit(x,y,1);
f = polyval(p,x);

% regression line values are calculated
ydiff = y-f;
xmax = max(x);
xmin = min(x);
%100 is arbitrary factor to give smooth plot of the fit
xx = (xmin:xmax/100:xmax); 
g = polyval(p,xx);

% the degrees of freedom are found
[r,c] = size(x);
N = r; % number of data pairs
nu = r-2; % degrees of freedom

% warning message when nu < 1
if nu<1
   disp(' ')
   disp('WARNING: degrees of freedom less than one >> invalid results !!!')
end

% various quantities are calculated
syx = sqrt(1/nu)*norm(ydiff,2);
xmean = mean(x);
ymean = mean(y);
ymdiff = y-ymean;
xmdiff = x-xmean;
sxx = sum((x-xmean).*(x-xmean));
a = norm(ydiff)^2;
b = norm(ymdiff)^2;
regco = sqrt(1-(a/b)); % linear regression coefficient
int = polyval(p,0); % intercept
slo = (polyval(p,xmax)-polyval(p,0))/xmax; % slope
tnup = tinv((1+.01*P)/2,nu); % student-t
preint = tnup*syx;
uinter = preint*sqrt((1/N)+(xmean^2/sxx)); % intercept uncertainty
uslope = preint*sqrt(1/sxx); % slope uncertainty
% JLO y2=x2*slo+int; % was used to scale y axis uncertainty in x-from-y plot

% compute the precision intervals for case 1
gu = g+preint;
gl = g-preint;

% compute the precision intervals for case 3
n = 1; % this sets the number of replicates of y used to compute y average
gamma = sqrt((1/n)+(1/N)+((xx-xmean).*(xx-xmean)/sxx));
grepu = g+preint*gamma;
grepl = g-preint*gamma;

% the data is plotted with y-error bars
subplot(2,1,1)
errorbar(x,y,ey,'ko');
xlabel('x'); ylabel('y');
axis([xmin xmax min(y) max(y)])
title(['P = ',num2str(P),'%; N = ',num2str(r),'; n = ',num2str(n),';S_{yx} = ',num2str(syx),';u_{int} = +-',num2str(uinter),'; u_{slo} = +-',num2str(uslope)])
hold on

% the linear regression line is plotted
plot(xx,g,'-k')
hold on

if caseno1==1
   plot(xx,gu,':r',xx,gl,':r')
   hold on
end

if caseno2==1
   plot(xx,grepu,':g',xx,grepl,':g');
   hold on
	% the uncertainty in x using the x-from-y estimate is determined
	% based upon the mean value of y and the precision intervals for case 3
   
   low=floor(min(y)); hi=ceil(max(y));
   for ypick=low:abs(hi-low)./xes:hi
   xforypick=(ypick-int)/slo;
	tol=0.00001; % convergence tolerance
	% find xg1 using Newtonian iteration
	xg1=xforypick; 
	xold1=xg1-abs(xmax-xmin);
	h=0.01;
	ct1=0;
	while abs(xg1-xold1)>tol
   	ct1=ct1+1;
   	xold1=xg1;
   	zeta1 = sqrt((1/n)+(1/N)+((xold1-xmean)^2/sxx));
   	zeta1h = sqrt((1/n)+(1/N)+((xold1-h-xmean)^2/sxx));
		F=xold1+preint*(zeta1/slo)-xforypick;
   	Fp=[(xold1+preint*(zeta1/slo)-xforypick)-...
   	    (xold1-h+preint*(zeta1h/slo)-xforypick)]/h;
		xg1=xold1-(F/Fp);
   	if ct1>100 break; end
	end
	% find xg3 using Newtonian iteration
	xg3=xforypick;
	xold3=xg3-abs(xmax-xmin);
	ct3=0;
	while abs(xg3-xold3)>tol
   	ct3=ct3+1;
   	xold3=xg3;
   	zeta3 = sqrt((1/n)+(1/N)+((xold3-xmean)^2/sxx));
		zeta3h = sqrt((1/n)+(1/N)+((xold3-h-xmean)^2/sxx));
		G=xold3-preint*(zeta3/slo)-xforypick;
		Gp=[(xold3-preint*(zeta3/slo)-xforypick)-...
   	    (xold3-h-preint*(zeta3h/slo)-xforypick)]/h;
		xg3=xold3-(G/Gp);
   	if ct3>100 break; end
	end
	delxlo=abs(xforypick-xg1);
	delxhi=abs(xg3-xforypick);
	delx=max(delxlo,delxhi);
   frac=100*delx/xforypick;
  	uncx=[uncx xforypick]; uncy=[uncy frac];
   end
   subplot(2,1,2)
	plot(uncx(:,2:end),uncy(:,2:end),'ok'); % JLO
	xlabel('x'); ylabel('% uncertainty in x');
    % JLO  axis([xmin xmax 0 y2]);
    axis([xmin 1.1*xmax 0 1.1*max(uncy(:,2:end))]); 
    title(['uncertainty in x-from-y estimate']);
end