clear
echo on
%Solution to problem 15
%We are asked to examine the production statistics
%of pecans, almonds, and walnuts.  The data from
%the problem has been stored in the file nuts.dat
%(actually, clip off the data and store it under
%this name in your directory).  Let's load this:
load 'nuts.dat'
nuts
pause
%We break this apart into two arrays:
year=nuts(:,1);
weight=nuts(:,2:4);
%The weight of the pecans is in the first column,
%that of almonds in the second, and walnuts in the
%third.
pause
%Let's get the means and the variances:
npt=length(year);
mean=sum(weight)/npt
sqr=sum(weight.*weight)/npt;
variance=npt/(npt-1)*(sqr-mean.^2)
%We can plot up this data:
errorbar([1:3],mean,variance.^.5,variance.^.5,'ro')
xlabel('type of nut')
ylabel('weight produced, kton')
axis([.5 3.5 0 600]);
hold on
plot(weight','grx')
hold off
pause
%Here the mean is given by a red circle, the spread
%due to the standard deviation by the errorbars
%and the data points by green x's.
%
pause
%Now for the covariance (actually, this bit of code gets
%the variances as well):
for i=1:3
 for j=1:3
  var(i,j)=npt/(npt-1)*(sum(weight(:,i).*weight(:,j))/npt-mean(i)*mean(j));
  echo off
 end
end
echo on
var
pause
%Note that this matrix is symmetric, as it should be, and
%that all of the covariances (the off-diagonal elements)
%are positive.  This is reasonable in that seasonal fluctuations
%lead to an increase or decrease in all nut production, at
%least on average.
%
pause
%Now let's compute the fraction pecans make of the total nuts
%produced.  In this part we only use the information on the
%mean and the variances and covariances:
fraction=mean(1)/sum(mean)
%That part was easy.  Now for the hard part.  We need to
%take the derivative of this function (fraction) with respect
%to the weight of pecans, almonds, and walnuts.  We can do this
%numerically:
eps=1e-8;
for i=1:3
 tmean=mean;
 tmean(i)=tmean(i)+eps;
 tfrac=tmean(1)/sum(tmean);
 deriv(i)=(tfrac-fraction)/eps;
 echo off
end
echo on
pause
%Now the derivatives are stored in the array deriv:
deriv
pause
%The variance can be captured by the simple formula:
fracvariance=deriv*var*deriv'
%which is exactly equivalent to the formula shown in class!
%
pause
%Now let's look at the last part.  We can calculate the
%fraction directly from the data:
%
frac=weight(:,1)./sum(weight')'
%
%This got a little tricky - I used the fact that sum sums over
%the columns of its arguments, so I transposed 'weight' before
%I summed over it, and then I transposed it back again.  This is
%equivalent to summing over each row.
pause
%Now to finish the problem:
exactfrac=sum(frac)/npt
exactvar=npt/(npt-1)*(frac'*frac/npt-exactfrac^2)
%which compares closely to the value calculated using the
%covariances.  It is not exact because the errors in the
%weights of the nuts are not particularly small relative to
%the means, and thus the approximate formula is in error.

%As a last note, we have only looked at population means and
%variances in this program, rather than the variance or error
%in the mean weights and mean ratio.  The problem as stated
%was a bit unclear as to which error it requested.  The two
%differ by just the square root of the number of points, e.g.
%just npt^.5   For this type of problem the population variance
%is probably more appropriate.  I would accept either answer,
%although you should always be careful in stating which error
%you are presenting.