clear
echo on
%Solution to Quiz # 4

%This script uses the attached function "radgrad.m" to analyze the data of Karnis, Goldsmith
%and Mason on concentration distributions in Poiseuille flow.  We have the measured array
%of bin counts:

n=[84 76 69 71];

%The calculated value of r*dc/dr is just:
ratio=radgrad(n)

%Which is a lot smaller in magnitude than the expected value of -0.135.  What is the error?

%We need the gradient with respect to n.  We do this numerically:

gradient=zeros(size(n));
for i=1:length(n)
 ntemp=n;
 ntemp(i)=ntemp(i)+1;
 gradient(i)=radgrad(ntemp)-ratio;
 echo off
end
echo on

%And that's it!  We adjusted each of the "n" values by one, sufficient
%to get an estimate of the gradient.

%We have the variance in n:
varn=diag(n);

%So the variance in the ratio is:
varratio=gradient*varn*gradient'

%Which yields the 95% confidence interval:
interval=[ratio-2*varratio.^.5,ratio+2*varratio.^.5]

%Which doesn't quite contain the expected value of -0.135, but is quite close...

%OK Question 2: There is no "right answer" to this one, but there are a few possibilities.
%
%First, the expected gradient ratio really isn't very far from the 95% confidence interval
%for these experiments.  There is a significant chance (around 1% to 2%) of them getting
%these results by random chance.  In fact, this sort of thing happens all the time: effects
%seem significant in one experiment, but are later shown not to be there after all, either
%because of randomness, or because of some missing factor.
%
%Second, this sort of measurement is extremely sensitive to operator error: The operator
%has to "decide" if a particle is really "sharply in focus".  This is even trickier for
%this particular experiment because the optical clarity is, in general, worse in the center
%of the tube rather than nearer the walls, because you are looking through more of the
%suspension to get to the central plane.  There may be some bias here which reduces the
%apparent concentration at the centerline.  Note that this is -not- taken into account in
%making the error estimate above, because it is systematic error rather than random.
%
%Third, although the size of the data set is much too small to be sure, the
%degree of variability in the data is actually much -less- than would be
%expected from Poisson statistics.  To see this, we may estimate the
%population standard deviation of the data from the "before and after" data
%for each bin.  We have the "before" data:
nb=[80 76 72 71];
%If we assume that there was -no- variation in the concentration in time,
%we get a standard deviation:
nmean=(nb+n)/2;
nvar=(nb-nmean).^2+(n-nmean).^2; %Note that we divide by 2-1=1 here...
npopstdev=mean(nvar)^.5
%which is quite a bit less than that expected from Poisson statistics:
nsdexpected=mean(nmean)^.5
%which is rather worrisome...
%
%The conclusion from this is that it is necessary to both do the statistics right, and to
%be very suspicious of the final results...which is also why it is important for results
%to be replicated by independent laboratories.
echo off