clear
echo on
%Solution to problem 4
%
%In this problem we wish to use the model
%and data given in class to predict the
%position of a lead weight 5 seconds after
%it is tossed up into the air.  The model
%is of the form:
%
%  h = x(1) + x(2)*t + x(3)*(-t^2/2)
%
%We first set up the linear regression problem:
pause
%We have the data:
t=[0;1;2;3;4];
h=[0.3;14.8;20.7;15.9;2.0];
[t,h]
pause
%This gives the matrix of modeling parameters A:
n=length(t)
a=[ones(n,1),t,-t.*t/2]
pause
%The solution vector is just:
x=a\h
%where the third modelling parameter is the acceleration
%due to gravity.
pause
%the position at t=5sec is just:
pvec=[1,5,-12.5];
pos=pvec*x
%and the velocity is just:
velvec=[0,1,-5];
vel=velvec*x
pause
%Now we have to calculate the error in these values.
%We shall do this using the normal equation formulation.
%We have that x = inv(A'*A)*A' * h, or, equivalently,
%x = K * h where K is appropriately defined.  Thus:
k=inv(a'*a)*a'
%which is a 3 x 5 matrix.  The first row gives the first
%coefficient, the second the second coefficient, etc.
pause
%We first have to calculate the variance in the measured
%positions h.  This is done from the residual:
m=length(x);
r=a*x-h;
varh=r'*r/(n-m);
stdh=varh^.5
%which is fairly small.  Remember that we have made the
%two key assumptions that the error in each measured
%value h is independent, and that they are all characterized
%by the same variance.  This means that the matrix of
%covariance of the measured values is just a single
%scalar variance times the identity matrix.
pause
%Next we calculate the matrix of covariance of the fitting
%parameters.  These can be calculated from the formula:
varx=k*k'*varh
pause
%We recognize that the position of the weight can be
%represented as a linear combination of fitted parameters:
%pos = pvec * x
%thus we have:
pvar=pvec*varx*pvec';
pstd=pvar^.5
%so the error in the predicted position is actually larger than
%the error in the measured values!  This is because the predicted
%position is dominated by the final measured value.  We can see
%this by looking at pvec*k:
pause
pvec*k
%If we had chosen an intermediate time (say 2.5 seconds) the error
%would have been much smaller.
pause
%The error in velocity can be calculated similarly:
velvar=velvec*varx*velvec';
velstd=velvar^.5
%which is also not all that small.
pause
%We can plot up the uncertainty interval of the correlation.  To
%do this we calculate the uncertainty of the predictions at the
%positions tp=[-1:.1:5].  Thus:
tp=[-1:.1:5]';
np=length(tp);
ap=[ones(np,1),tp,-tp.*tp/2];
ppos=ap*x;
varppos=diag(ap*varx*ap');
%where we are only interested in the variance of each
%position rather than the covariances.
pause
%Now for the plotting part:
%It turns out that the error is very small.  To be able to see
%the difference between the upper and lower bounds of the
%predicted positions, we shall multiply the error by 10.  Thus:
figure(1)
plot(t,h,'o')
xlabel('time')
ylabel('height')
hold on
plot(tp,ppos,'g')
plot(tp,ppos+10*varppos.^.5,'--r')
plot(tp,ppos-10*varppos.^.5,'--r')
text(1,-20,'Error bounds exaggerated by 10x')
hold off
pause
%Note that the magnitude of the error grows as we move away
%from the region where the measurements are made.
echo off