clear
echo on
%Solution to Quiz #3
%We have the average count rates:

s=35
s0=40
b=32

%We also have the decay rate:
lambda=log(2)/5730

%OK, now for question 1:
age=-log((s-b)/(s0-b))/lambda

%Question 2 is a little harder.  First we need the variances in count rates:
t=400 %The length of time we are counting...
vars=s/t
vars0=s0/t
varb=b/t

%and we need the gradient.  It is easiest to do this numerically, actually.
%We use a forward difference algorithm:
ep=0.001
grads=(-log((s+ep-b)/(s0-b))/lambda-age)/ep;
grads0=(-log((s-b)/(s0+ep-b))/lambda-age)/ep;
gradb=(-log((s-b-ep)/(s0-b-ep))/lambda-age)/ep;

%So we get the gradient vector:
gradage=[grads,grads0,gradb]

%The three measurements are independent (no covariance), so the variance in 
%the age is:
varage=gradage*diag([vars,vars0,varb])*gradage'

%and the standard deviation:
stdage=varage.^.5

%Now for question 3:
%We calculate the confidence interval for the ratio (S-B)/(S0-B):
ratio=(s-b)/(s0-b)
%we have the gradient:
gradratios=((s+ep-b)/(s0-b)-ratio)/ep;
gradratios0=((s-b)/(s0+ep-b)-ratio)/ep;
gradratiob=((s-b-ep)/(s0-b-ep)-ratio)/ep;
gradratio=[gradratios,gradratios0,gradratiob]

sdratio=(gradratio*diag([vars,vars0,varb])*gradratio')^.5

%The 95% confidence interval is +/- 2 sigma (for a normal distribution), thus
ratiointerval=[ratio+2*sdratio,ratio-2*sdratio]

%which yields the age interval:
ageinterval=-log(ratiointerval)/lambda

%or, subtracting off the mean age:
deviation=ageinterval-age

%You can see that we have a pretty large standard deviation!  Also, it isn't
%uniformly distributed around the mean due to the non-linearity.

%Finally, for Question 4:
%The gradient is unchanged by the count time - only the variances of the 
%measured count rates.  The standard deviation needs to be reduced by a factor
%of about:
factor=(deviation(2)-deviation(1))/400

%We can get this if we increase the time by this factor squared!  Thus:
newtime=t*factor^2

%or about a month...  Needless to say, I didn't get very accurate age
%estimates...
echo off