% PS02_3.m is the m-file answer for problem set 2 problem 3. It was created 
% by T. Kenna 980924 1222
% Modif'd 2002:09:25 DMG to fix Tim's errors in propagation of variance and his
%                    backward sense of He:Ne ratios

clear
close all
t=10;			% The display time for figures

% clc
ans1=['Let''s start by loading the data files He.dat and Ne.dat.  I know'];
ans2=['you weren''t asked to plot, but its not a bad idea to see what your'];
ans3=['data looks like. Again, any figures I create will only live for 10'];
ans4=['seconds, if you need more time edit the program and change the value'];
ans5=['of t. <return>'];
ans6=[' '];
disp(ans1),disp(ans2),disp(ans3),disp(ans4),disp(ans5),disp(ans6)
pause

load He.dat             % Load the data
load Ne.dat		% Load the data
figure(1);clf 		% Select figure 1 and clear any previous stuff 	
plot(He,Ne,'bo') 		% Plot the data
title('PS02-3 He vs Ne Plot')
xlabel('He')
ylabel('Ne')
pause(t)		% View the plot for a limited time
delete(figure(1))	% Kill it and return to matlab

% clc			% clear the text
ans1=['Okay, now that we have seen our data, let''s do the the two type I'];
ans2=['regressions using the commands [a,sa,cov,r] = linfit(He,Ne,0) and'];
ans3=['[a,sa,cov,r] = linfit(Ne,He,0).  Remember, the first variable is the'];
ans4=['independent variable, the second is the dependent variable, and the'];
ans5=['third variable is the uncertainty in the dependent variable.  If you'];
ans6=['type help linfit, you will see that a single number in this spot'];
ans7=['rather than a vector the size of your dependent variable will result'];
ans8=['in an unweighted fit.  Since we can distinguish between our two'];
ans9=['diffusion models by the He:Ne ratio, we need only be concerned with'];
ans10=['the slopes and their uncertainty (ie a(2) and sa(2) for each'];
ans11=['regression). <return>'];
disp(ans1),disp(ans2),disp(ans3),disp(ans4),disp(ans5),disp(ans6),disp(ans7)
disp(ans8),disp(ans9),disp(ans10),disp(ans11)
pause

% clc
[a,sa,cov,rhn] = linfit(Ne,He,0);		% The first type I He on Ne
m1=a(2);				% Save the slope
sm1=sa(2);				% Save the uncertainty
He_Ne_Type_1=[m1 sm1]			% Package the results
[a,sa,cov,rnh] = linfit(He,Ne,0);		% The second type I Ne on He
m2=a(2);				% Save the slope
sm2=sa(2);				% Save the uncertainty
Ne_He_Type_1=[m2 sm2]			% Package the results
ans1=['Mmm, both regressions give a slope of ~0.89; ideally, they should be'];
ans2=['the inverse of one another.  The optimal fit is obtained by taking'];
ans3=['the geometric mean of the two with the command sqrt(m1/m2). In'];
ans4=['addition, we can propagate the uncertainty in the geometric mean'];
ans5=['using the command sqrt(m1/m2)*sqrt(sm1^2/m1 + sm2^2/(4*m2^2)) <return>'];
%% sqrt((sm1/m1+sm2/m2)^2)*sqrt(m1/m2). 
disp(ans1),disp(ans2),disp(ans3),disp(ans4),disp(ans5)
pause

% clc
Geometric_mean=[sqrt(m1/m2) sqrt(m1/m2)*sqrt(sm1^2/m1 + sm2^2/(4*m2^2))]
ans1=['But keep in mind that the error estimate is probably inaccurate'];
ans2=['because we have NO WAY of estimating the covariance term. Certainly'];
ans3=['m1 and m2 covary, so this is just a ... guesstimate. <return>'];
ans4=[' '];
disp(ans1),disp(ans2),disp(ans3),disp(ans4)
pause

ans1=['Next, let''s do the regression as fully type II (the more correct'];
ans2=['approach) with the command [m,b,r,sm,sb,xbar,ybar] = lsqfitma(He,Ne)'];
ans3=['and compare the slopes and uncertainties.  <return>'];
disp(ans1),disp(ans2),disp(ans3)
pause

% clc
[m,b,r,sm,sb,xbar,ybar] = lsqfitma(Ne,He);	% The type II regression
Slope_comparison=[Geometric_mean;m sm]		% The comparison
ans1=['Both approaches yield essentially the same value. This slope (which'];
ans2=['is the ratio of He:Ne in bubble injection/gas exchange, is within'];
ans3=['errors of the second model (which predicts a slope of 1.027) but is'];
ans4=['not consistent with the first model (which predicts a slope of'];
ans5=['0.956). Our estimated slope from the type II regression is about 1.6'];
ans6=['standard deviations from the latter, suggesting that we can reject'];
ans7=['the first model at about the 97% confidence level. Whereas the our'];
ans8=['Type-II slope is only about 1.1 standard deviations from the second'];
ans9=['model. We can reject the second model with only 93% confidence level.'];
ans10=['<return>'];
disp(ans1),disp(ans2),disp(ans3),disp(ans4),disp(ans5),disp(ans6),disp(ans7),disp(ans8),disp(ans9),disp(ans10)
pause

% clc
ans1=['End of problem 3'];
disp(ans1)