clear
echo on
%Solution to Problem 6
%
%Save the data file heat.dat in the appropriate
%directory before running this program!
%
%In this problem we were asked to analyze
%some heat transfer data.  We load it in:
load heat.dat
re=heat(:,1);
pr=heat(:,2);
nu=heat(:,3);
pause
%We wish to fit the data to a power law
%model of the form:

%Nu = c * Re^a * Pr^b

%To use linear regression we must linearize
%the problem.  We do this by taking the log.
%Thus we have the regression problem:
n=length(nu);
m=3;
a=[log(re),log(pr),ones(n,1)];
b=log(nu);
%which has the solution:
x=a\b
%where x(1) is a, x(2) is b, and x(3) is log(c).
pause
%We may also calculate the matrix of covariance:
k=inv(a'*a)*a';
r=b-a*x;
varb=r'*r/(n-m)
varx=varb*k*k'
pause
%We thus have the values for x and the errors:
[x,diag(varx).^0.5]
%Note that while we can calculate the uncertainties in
%each of these fitting parameters, they are NOT
%independent variables, as all arise from the same
%data set!  This is reflected in the off-diagonal
%elements of the matrix of covariance!
pause
%Now we do it again using the bootstrap method.
%We shall take p samples of n data points from
%our replicated data set.  We do this using the
%random number generator provided by matlab:
p=500;
for i=1:p
 index=ceil(rand(n,1)*n);
 bboot=b(index);
 aboot=a(index,:);
 xboot(:,i)=aboot\bboot;
 echo off
end
echo on
pause
%We now have an array of solutions for our random
%samples of our replicated data set.  Let's look at
%the average of x:
avexboot=(sum(xboot')')/p
%which is close to the value which we had before:
ratio=avexboot./x
pause
%And then the matrix of covariance:
for i=1:m
 for j=1:m
  varxboot(i,j)=(xboot(i,:)*xboot(j,:)'/p- ...
  avexboot(i)*avexboot(j))*p/(p-1);
  echo off
 end
end
echo on
varxboot
pause
%Finally, we compare this to varx:
ratiovar=varxboot./varx
%which matches pretty closely.  Remember that
%the ratio of the standard deviations are the 
%square roots of these values.  Note that the
%value of p which we have chosen is much larger
%than n.  We need to have a large value for p
%because it is much harder to estimate the value
%of the variance than it is to estimate the mean
%(you need alot more points).  Remember that
%because the bootstrap uses a random sample, it
%(the means and variances) will change every
%time you run it!
pause
%It is also useful to do a little graphics.  This is
%tricky, because we have two independent variables:
%Re and Pr.  In fact, however, there is very little
%change in Pr, and the dependence on Pr is quite
%weak (since we aren't changing the fluid, it only
%changes due to different temperatures).  This is
%also reflected in the fairly large error in the
%calculated exponent.  We can look at this graphically.
pause
%Let's plot up the data first by simply ignoring
%the Pr dependence.  We give the "line" by taking
%the correlation for the average Pr of the exp'ts:
avepr=mean(pr)
figure(1)
loglog(re,nu,'o')
hold on
plot(re,exp(x(1)*log(re)+x(2)*log(avepr)+x(3)),'r')
hold off
xlabel('Re')
ylabel('Nu')
legend('uncorrected data','model')
pause
%We can also "correct" for the calculated pr dependence
%by dividing the Pr dependence out of the data:
nuadj=nu./(pr/avepr).^x(2);
figure(2)
loglog(re,nu,'o',re,nuadj,'*')
hold on
plot(re,exp(x(1)*log(re)+x(2)*log(avepr)+x(3)),'r')
hold off
xlabel('Re')
ylabel('Nu')
legend('uncorrected data','adjusted data','model')
title(['Heat transfer data adjusted to Pr = ',num2str(avepr)])
%As you can see, it reduces the scatter a little bit, but
%not very much.
pause
%We can go the other way, extracting the Re dependence to
%look at how it varies with Pr - but we expect this to be
%-really- rough, as there just isn't much Pr variation.
avere=mean(re)
nuadj2=nu./(re/avere).^x(1);
figure(3)
loglog(pr,nuadj2,'o')
hold on
plot(pr,exp(x(1)*log(avere)+x(2)*log(pr)+x(3)),'r')
hold off
xlabel('Pr')
ylabel('Nu')
legend('adjusted data','model')
title(['Heat transfer data adjusted to Re = ',num2str(avere)])
%The data is pretty noisy this way (which is why you get
%a large error on the exponent), but you can pick off the
%positive dependence on Pr.
echo off