%% File census_model.m
%
% This uses the 11x2 census data matrix 
%
census = [1900. 77584.E3;
          1910. 92407.E3;
          1920. 106461.E3;
          1930. 123077.E3;
          1940. 132122.E3;
          1950. 152271.E3;
          1960. 180671.E3;
          1970. 205052.E3;
          1980. 227224.E3;
          1990. 249438.E3;
          2000. 282171.E3];
%
% to do the following:
%
% (1) Prepares three 11x11 matrices A, B, C corresponding to 
%     10th degree polynomial models in t, (t-1940), and (t-1940)/40;
%     prints out condition numbers of the LS problems associated
%     with polynomials of degrees 0 through 10 for each model.
%
% (2) Graphs the LS polynomial models of degrees 2:2:10 in (t-1950)/50 
%     for 1890 <= t <= 2010, superimposing the data in each case.  
%
% (3) Shows the effect on the 10th-degree polynomial of perturbing the
%     data by -1.e6 in the direction of the singular vector associated 
%     smallest singular value. 
%
% Note: These census data were obtained from 
% https://fusiontables.google.com/DataSource?dsrcid=225439#rows:id=1
%

b=census(:,2); 
z=census(:,1); 

% Plot the data. 
figure(1);plot(z,b,'*r');axis([1890 2010 0 3.e08]);
title({'Census Data, 1900 - 2000'; 'Strike any key to continue.'}); drawnow;
pause;

% Set up the matrices A, B, C for the LS problems. 
A=[ones(11,1),z,z.^2,z.^3,z.^4,z.^5,z.^6,z.^7,z.^8,z.^9,z.^10];
y=z-1950;
B=[ones(11,1),y,y.^2,y.^3,y.^4,y.^5,y.^6,y.^7,y.^8,y.^9,y.^10];
w=y/50;
C=[ones(11,1),w,w.^2,w.^3,w.^4,w.^5,w.^6,w.^7,w.^8,w.^9,w.^10];
% Print out the condition numbers of the least-squares problems. 
fprintf('\n\n Condition numbers of the 11xk LS problems:\n')
fprintf('  k        t           (t-1950)      (t-1950)/50\n')
for k=1:11
  fprintf('  %d   %e   %e   %e \n', k, cond(A(:,1:k)), ...
     cond(B(:,1:k)), cond(C(:,1:k))); 
end
fprintf('\n Strike any key to continue.\n');pause;

x = zeros(11,1);
t = [1890:1:2010]';
D = zeros(size(t,1),11);
D(:,1) = ones(size(t,1),1);
for k = 2:11
  D(:,k) = D(:,k-1).*((t-1950)/50);
end
for kk = 1:5
  k = 2*kk+1;
  x(1:k,1) = C(:,1:k)\b;
  u = D(:,1:k)*x(1:k,1);
  figure(kk);
%  subplot(4,1,kk);
  plot(t,u);hold on;plot(z,b,'*r');axis([1890 2010 0 3.e08]);drawnow;
  title({['Polynomial Degree ', num2str(k-1)];'Strike any key to continue.'});hold off;
  pause;
end
% Uncomment the lines below to see effect of perturbing along singular 
% vector corresponding to smallest singular value. 
[U,S,V] = svd(C);delta = -1.e6; b = b + delta*U(:,11);
x = C\b;u = D*x;
figure(gcf);
plot(t,u);hold on;plot(z,b,'*r');axis([1890 2010 0 3.e08]);drawnow;
title(['Polynomial Degree ', num2str(k-1), ' with Perturbed Data.']);
hold off;
