home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / CENSUS.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.3 KB  |  129 lines
  1. %CENSUS    Try to predict the US population in the year 2000.
  2.  
  3. %    C. Moler & B. Jones, 10/27/92.
  4. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  5.  
  6. fig = clf;
  7. clc
  8. echo on
  9.  
  10. % This example is older than MATLAB.  It started as an exercise in
  11. % "Computer Methods for Mathematical Computations", by Forsythe,
  12. % Malcolm and Moler, published by Prentice-Hall in 1977.
  13. %
  14. % Now, MATLAB and Handle Graphics make it much easier to vary the
  15. % parameters and see the results, but the underlying mathematical
  16. % principles are unchanged:
  17.  
  18. %    Using polynomials of even modest degree to predict
  19. %    the future by extrapolating data is a risky business.
  20.  
  21. pause  % Press any key to continue after pauses.
  22.  
  23. clc
  24.  
  25. % Here is the US Census data from 1900 to 1990.
  26.  
  27. % Time interval
  28. t = (1900:10:1990)';
  29.  
  30. % Population
  31. p = [75.995 91.972 105.711 123.203 131.669 ...
  32.      150.697 179.323 203.212 226.505 249.633]';
  33.  
  34. figure(fig);
  35. plot(t,p,'go');
  36. axis([1900 2010 0 400]);
  37. title('Population of the U.S. 1900-1990');
  38. ylabel('Millions');
  39.  
  40. pause % Press any key.
  41.  
  42. clc
  43.  
  44. % What is your guess for the population in the year 2000?
  45. p
  46.  
  47. pause % Press any key.
  48.  
  49. clear A
  50. clc
  51.  
  52. % Let's fit the data with a polynomial in t and use it to
  53. % extrapolate to t = 2000.  The coefficients in the polynomial
  54. % are obtained by solving a linear system of equations involving
  55. % a 10-by-10 Vandermonde matrix, whose elements are powers of
  56. % scaled time, A(i,j) = s(i)^(n-j);
  57.  
  58. n = length(t);
  59. s = (t-1900)/10;
  60. A(:,n) = ones(n,1);
  61. for j = n-1:-1:1
  62.    A(:,j) = s .* A(:,j+1);
  63. end
  64.  
  65. pause  % Press any key.
  66.  
  67. clc
  68.  
  69. % The coefficients c for a polynomial of degree d that fits the
  70. % data p are obtained by solving a linear system of equations
  71. % involving the last d+1 columns of the Vandermonde matrix:
  72. %    
  73. %        A(:,n-d:n)*c ~= p
  74. %
  75. % If d is less than 9, there are more equations than unknowns and
  76. % a least squares solution is appropriate.  If d is equal to 9, the
  77. % equations can be solved exactly and the polynomial actually interpolates
  78. % the data.  In either case, the system is solved with MATLAB's backslash
  79. % operator.  Here are the coefficients for the cubic fit.
  80.  
  81. c = A(:,n-3:n)\p
  82. pause  % Press any key.
  83.  
  84. clc
  85.  
  86. % Now we evaluate the polynomial at every year from 1900 to 2010,
  87. % plot the results, and save the "handles" of the plots for later use.
  88.  
  89. v = (1900:2010)';
  90. x = (v-1900)/10;
  91. y = polyval(c,x);
  92. z = polyval(c,10);
  93.  
  94. figure(fig);
  95. hold on
  96. yhandle = plot(v,y,'m-');
  97. zhandle = plot(2000,z,'yx');
  98. ztext = text(2000,z-10,num2str(z));
  99. title('Polynomial fit, degree = 3')
  100. hold off
  101.  
  102. pause  % Press any key.
  103.  
  104. clc
  105.  
  106. % Finally, we establish a vector of Handle Graphics buttons
  107. % that allow the degree of the polynomial to be varied.
  108. % The details are in the M-files buttonv.m and censusex.m.
  109.  
  110. d = 3;
  111. buttons = buttonv([.18 .8], 0:9, d, 'censusex', ...
  112.    'spline', 'exp', 'bounds', 'done');
  113.  
  114. % We also provide interpolation by a cubic spline, a least squares
  115. % fit by an exponential in t and optional confidence bounds.  The
  116. % confidence bounds account for statistical errors in the data,
  117. % but assume the polynomial model is correct.
  118.  
  119. echo off
  120.  
  121. % Establish handle for bounds, but don't plot anything yet.
  122. bounds = 0;
  123. figure(fig);
  124. hold on
  125. boff = [NaN*y; NaN; NaN*y];
  126. bhandle = plot([v; NaN; v], boff,'c');
  127. offscale = 0;
  128. hold off
  129.