% Sample solution for final exam problem in MAT 342
% Linear Algebra, Summer 2000. Matthias Kawski
% http://math.la.asu.edu/~kawski
%
% Open the source file in the editor to read detailed
% comments using:   > edit logpopsol
 
%Read in the supplied data, display and plot them
[t,y]=logpop;
[t,y]
plotlog=figure
set(plotlog,'Position',[50 300 300 300])
plot(t,y,'ro')
axis([-2,20,1,6])
text(2,1.3,'Time: Decades since 1800')
text(-1,2,'Log of population (millions)','Rotation',90)
disp('Strike any key to continue')
pause
%
% For a quadratic fit create the matrix of function
% values of the basis functions {1,t,t^2} at given 
% times t(i).
%
A=[ones(size(t)),t,t.*t]
disp('Strike any key to continue')
pause
%
% Solve the normal equations (since the problem size 
% is small, there is no need for QR-factorizations...)
%
c=(A'*A)\(A'*y)
%
% Compare with the MATLAB supplied default least squares 
% solution
%
cc=A\y
%
% The solution: the components of c are the coordinates  
% of the desired best fitting parabola with respect to 
% the basis {1,t,t^2}.
disp('Strike any key to continue')
pause
%
yy=A*c;
disp(' ')
disp('     time      data       fit')
[1800+10*t,y,yy]
%
% Overlay the graph for an intuitive check of goodness 
%
hold on
figure(plotlog)
plot(t,yy,'b-')
text(1,5.5,'Quadratic model')
text(1,5.1,'log(pop)=a+b*y+c*t^2') 
text(9,3.5,'Best choice') 
text(9,3.2,'of parameters:') 
text(9,2.5,strcat('a=',num2str(c(1)))) 
text(9,2.2,strcat('b=',num2str(c(2)))) 
text(9,1.9,strcat('c=',num2str(c(3)))) 
%
disp('Strike any key to continue')
pause
%
% It is instructive to compare this solution with the
% least squares fit of the original population data, one
% i.e. a quadratic with the exponetial of a quadratic
%
P=exp(y);
b=A\P;
PP=A*b;
YY=exp(yy);
[t,P,YY,PP]
plotpop=figure
set(plotpop,'Position',[150 300 400 300])
hold on
plot(1800+10*t,P,'k*')
plot(1800+10*t,YY,'b-')
plot(1800+10*t,PP,'r:')
text(1760,220,'Two different models for the US population')
text(1760,190,'stars: measured pop data 1790 - 1990')
text(1760,170,'blue: pop=exp(a+b*t+c*t^2)')%
text(1760,150,'red:    pop=   a+b*t+c*t^2')
disp('Strike any key to continue')
pause
%
moret=[19:30]';
B=[ones(size(moret)),moret,moret.*moret];
plotextra=figure
set(plotextra,'Position',[200 300 500 300])
hold on
plot(1800+10*t,P,'k*')
plot(1800+10*t,YY,'b-')
plot(1800+10*t,PP,'r:')
plot(1800+10*moret,exp(B*c),'c-')
plot(1800+10*moret,B*b,'m:')
text(1790,580,'Extrapolation of US population using different models')
text(1790,500,'stars:= measured data 1790 - 1990')
text(1790,450,'cyan:       pop=exp(a+b*t+c*t^2)')
text(1790,400,'magenta: pop =   a+b*t+c*t^2 ')
disp('Strike any key to delete all figures')
pause
delete(plotlog)
delete(plotpop)
delete(plotextra)