clear
echo on
%Solution to problem 5
%
%In this problem we are asked to do a linear
%regression fit of a polynomial to the diffusivity
%data measured for dilute suspensions undergoing
%shear.  The standard deviation of the measurements
%is given, so we shall use a weighted linear 
%regression procedure.  First we put in the data:
c=[.01,.025,.05,.075,.1,.15]';
diff=[2.37e-4,5.63e-4,1.4153e-3,2.27e-3,4.055e-3,9.965e-3]';
sigma=[2.2e-5,6.491228e-5, 1.617643e-4,3.98208e-4,7.605929e-4,1.775587e-3]';
pause
%Let's plot this data up:
figure(1)
errorbar (c,diff,sigma)
xlabel('concentration')
ylabel('diffusivity')
%As you can see, it certainly goes to zero as the
%concentration goes to zero.
pause
%The asymptotic behavior can be visualized better
%by plotting diff/c rather than just diff:
figure(2)
errorbar(c,diff./c,sigma./c,sigma./c,'o')
axis([0 .16 0 .1])
xlabel('concentration')
ylabel('diffusivity / concentration')
%Here there is clearly a non-zero asymptotic value
%of diff/c as c goes to zero.  This is what we are
%trying to capture.
pause
%OK, now for the weighted linear regression.  The
%unweighted linear regression problem is of the 
%form:
%
%   min || A x - b ||
%
%In the weighted regression problem we weight each
%data point by the inverse of its variance when forming
%the sum of the squares of the deviation.  Let's do this:
pause
%First we set up the matrix of modeling functions.  We
%have a linear term, a quadratic term, and a cubic term:
a=[c,c.^2,c.^3]
pause
%This is the weighting function:
vardiff=diag(sigma.^2);
weight=inv(vardiff)
pause
%We now define a matrix kw:
kw=inv(a'*weight*a)*a'*weight
%which can be used the conventional way:
x=kw*diff
%which has a lead coefficient different from zero.  Now
%for the error:
pause
varx=kw*vardiff*kw'
xerror=diag(varx).^.5;
pause
%So the coefficients and errors are:
[x,xerror]
%and in particular,
x(1)/xerror(1)
%So the order c coefficient is statistically different from
%zero (about 8 standard deviations).
pause
%Let's finish off with a bit of graphics.  Let's plot up
%the model predictions together with its uncertainty.
cp=[.001:.001:.15]';
ap=[cp,cp.^2,cp.^3];
dp=ap*x;
perror=diag(ap*varx*ap').^.5;
figure(2)
hold on
plot(cp,dp./cp,'r')
plot(cp,(dp+perror)./cp,'g')
plot(cp,(dp-perror)./cp,'g')
hold off
axis([0 .16 0 .1]);
%where we have plotted the 1-sigma confidence interval
%of the model function, based on the given standard deviations.

echo off