clear
reset(0)
echo on
format compact
%This program shows how you can use Matlab
%to do a simple error calculation.  The example
%is determining the acceleration due to gravity
%by measuring the height of a ball as a function
%of time which is initially thrown upward with 
%some velocity.  Thus:
pause
%First we set up the array of times:
t=[0:.05:2.5]';
%Note that we were careful to make this
%a column vector!

%Next we generate some artificial 'data'
%with random noise:
b=(-0.5*t.^2+1.25*t)*9.8+randn(size(t));

%and create the matrix of modelling functions:
a=[ones(size(t)),t,t.^2];
pause
%And this is the matrix k used in the solution:
k=inv(a'*a)*a';

%Here are the best fit modelling parameters:
x=k*b
pause
%Now we estimate the variance in b:
[n m]=size(a);
varb=norm(a*x-b)^2/(n-m);

%This yields the standard deviation in b:
stddev=varb^.5
%which is close to the expected value of unity
%for the randn routine.
pause
%We use the variance in b to calculate the variance
%in x:
varx=varb*k*k'
pause
%We want to generate a smooth curve for our
%best fit curve:
tp=[0:.1:2.5]';
ap=[ones(size(tp)),tp,tp.^2];
bp=ap*k*b;

%And we plot this up with the data:
figure(1)
plot(t,b,'o',tp,bp)
set(gca,'FontSize',18)
xlabel('time (sec)')
ylabel('height (m)')
drawnow
hold on
pause
%Now we calculate the variance in the model
%predictions:
varbp=varb*ap*k*k'*ap';

%And we extract the square root of the diagonal
%to get the standard deviation:
sigbp=diag(varbp).^.5;

%We add the 95% confidence interval to the plot:
plot(tp,bp+2*sigbp,'g--',tp,bp-2*sigbp,'g--')
hold off
pause
%And we finish by outputting the best fit parameter
%values and their standard deviations:
[x,diag(varx).^.5]

%The first coefficient is just half the acceleration
%due to gravity, thus the estimated gravitational
%acceleration, with random error is:
gravity = [-x(3),varx(3,3)^.5]*2

%which is reasonably close to the expected value of
%9.8 m/s^2.
pause
%Note that the data mostly lie outside the dashed lines.
%This is because the confidence interval is that of the
%-model-, specifically assuming that the error in the
%data is random.  Because we have alot of data points,
%the uncertainty in the model (under this assumption)
%is small.  If we had an infinite number of data points
%the calculated uncertainty would be -zero-, which is
%wrong!  In general, for alot of data points the error
%is dominated by systematic error, which is not captured
%by the assumptions.  We'll look at how to deal with this
%problem in the last lecture.
pause