function censusapp2024(arg)
%CENSUSAPP2024 Predict the US population beyond 2020.
% This example began as an exercise in "Computer Methods for
% Mathematical Computations", by Forsythe, Malcolm and Moler,
% published by Prentice-Hall in 1977. The software was written
% in Fortran, the predominate technical programming language at
% the time. Since 1970, the data has been updated every ten years.
% Today, MATLAB makes it easier to vary parameters and visualize
% results, but the underlying mathematical principles are unchanged:
%
% Using polynomials of even modest degree to predict
% the future by extrapolating data is a risky business.
%
% The data are from the decennial census of the United States for the
% years 1900 to 2020. The task is to extrapolate beyond 2020.
%
% In addition to polynomials of various degrees, you can choose
% interpolation by a cubic spline or two shape-preserving Hermite
% cubics, and least squares fits by an exponential, a logistic curve
% and an exponential of an exponential.
%
% spline C2, "not-a-knot".
% pchip C1, shape-preserving.
% makima C1, modified Akima.
%
% exponential exp(b*t+c)
% logistic a./(1+exp(-b*(t-c)))
% gompertz a*exp(-b*exp(-c*t))
%
% Error estimates account for errors in the data, but not in the model.
% R^2 is the proportion of variation in data predicted by the model.
%
% Usage:
% censusapp2024 Figure width 800.
% censusapp2024(figsize) Figure width 320 + 480*figsize.
% censusapp2024(method) Callback
%
% Date: 19-May-2024.
% Copyright 2014-2024 Cleve Moler
% Census data for 1900 to 2020.
t = (1900:10:2020)';
p = [ 75.995 91.972 105.711 123.203 131.669 150.697 ...
179.323 203.212 226.505 249.633 281.422 308.746 ...
331.449]';
x = (1890:2040)'; % Evaluation years
w = 2030; % Extrapolation target
z = 343; % Guess w population
d = 0; % Polynomial degree
r2 = 0; % R^2
level = 0; % Noise level
models = {'census','polynomial','pchip','spline','makima', ...
'exponential','logistic','gompertz'}';
if nargin == 0 || nargin == 1 && isscalar(arg)
if nargin == 0
figsize = 1;
else
figsize = arg;
end
model = 'census';
h = init_graphics;
noise = rand*randn(size(p));
set(gcf,'UserData',{p,h,noise,figsize});
else
model = arg;
fud = get(gcf,'UserData');
[p,h,noise,figsize] = deal(fud{:});
level = h.noise.UserData;
p = p + level*noise;
h.plot(1).YData = p;
end
% Update plot with new model
w = h.predict.UserData;
x = (1890:max(2040,w+30))';
switch model
case 'census'
y = NaN*x;
z = 324.790 + 2.64*(w-2017) + 0.0042*(w-2017)^2;
case 'polynomial'
[y,z,r2] = poly(h,t,p,x,w);
case {'pchip','spline','makima'}
[y,z] = hermite(model,t,p,x,w);
r2 = 1; % Exact fit
case {'exponential','logistic','gompertz'}
[y,z,r2] = expo(model,t,p,x,w);
end
update_graphics(model,h,x,y,w,z,r2)
% ------------------------------------------------
function h = init_graphics
figw = 320 + 480*figsize;
figh = 3/4*figw;
fontsize = 5 + 9*figsize;
fig = gcf;
fig.Position = [300 150 figw figh];
fig.Name = 'censusapp2024';
fig.MenuBar = 'none';
fig.NumberTitle = 'off';
shg
clf
h.plot = plot(t,p,'bo', x,0*x,'k-', ...
w,0,'.', [x;NaN;x],[x;NaN;x],'m:');
lw = 1+.06*(fontsize-5);
h.plot(1).LineWidth = lw;
h.plot(1).MarkerSize = fontsize-2;
h.plot(2).LineWidth = lw;
h.plot(4).LineWidth = lw;
h.plot(3).MarkerSize = 3*fontsize;
h.plot(3).Color = [0 2 0]/3; % darkgreen
h.plot(4).Color = [2 0 2]/3; % darkmagenta;
h.title = ['Predict U.S. Population in ' int2str(w)];
axis([1890 2060 0 450]);
x0 = .18; % Upper left corner of gui
y0 = .82;
dx = .08;
dy = .07;
dw = .05;
h.model = uicontrol('style','popup','string',models, ...
'units','normal','position',[x0 y0-3*dy .24 dw], ...
'fontsize',fontsize,'fontweight','bold','callback',@model_cb);
h.extrap = text(w,z+20,'predict','horiz','right', ...
'color',[0 2 0]/3,'fontsize',fontsize,'fontweight','bold');
h.predict = uibut('push',sprintf('predict = %d',w), ...
[x0+dw y0 .18 dw], [],[],w);
h.noise = uibut('push',sprintf('noise = %d',level), ...
[x0+dw y0-dy .18 dw], [],[],level);
h.deg = uibut('push',sprintf('degree = %d',d), ...
[x0+dw y0-2*dy .18 dw], [],[],d);
uibut('push','-', [x0 y0 dw dw], @predict_cb,[],[]);
uibut('push','+', [x0+3*dx y0 dw dw], @predict_cb,[],[]);
uibut('push','-', [x0 y0-dy dw dw], @noise_cb,[],[]);
uibut('push','+', [x0+3*dx y0-dy dw dw], @noise_cb,[],[]);
h.minus = uibut('push','-', [x0 y0-2*dy dw dw], @degree_cb,[],[]);
h.plus = uibut('push','+', [x0+3*dx y0-2*dy dw dw], @degree_cb,[],[]);
h.r2txt = uibut('push','R^2',[x0+.06 y0-4*dy .18 dw],[],[],1);
h.r2 = uibut('toggle','R^2',[x0 y0-4*dy dw dw],@r2_cb,1,[]);
h.err = uibut('toggle','errs',[x0 y0-5*dy dw .05],@errs_cb,0,[]);
if figsize >= 0.5
uibut('push','info',[.91 .90 .08 .06],@info_cb,[],[]);
uibut('push','help',[.91 .82 .08 .06],@help_cb,[],[]);
uibut('push','restart',[.91 .74 .08 .06],@restart_cb,[],[]);
ylabel('Millions')
end
end
function update_graphics(model,h,x,y,w,z,r2)
% Update plot
h.plot(2).XData = x;
h.plot(2).YData = y;
h.plot(3).Visible = 'on';
h.plot(3).XData = w;
h.plot(3).YData = z;
h.title = title(['Predict U.S. Population in ' int2str(w)]);
h.title.FontSize = h.predict.FontSize+2;
ax = axis;
h.extrap.Position = [w,min(max(z+20,30),ax(4)-30)];
h.extrap.Visible = 'on';
h.extrap.String = sprintf('%6.1f',z);
% Update controls
switch model
case 'census'
h.extrap.String = 'predict';
end
% Display error estimates if requested
if h.err.Value == 1
errest = errorestimates(model,t,p,x,y);
h.plot(4).Visible = 'on';
h.plot(4).XData = [x;NaN;x];
h.plot(4).YData = errest;
else
h.plot(4).Visible = 'off';
end
if h.r2.Value == 1
h.r2txt.String = sprintf('%8.6f',r2);
else
h.r2txt.Visible = 'off';
end
end
function h = uibut(style,string,position,callback,value,userdata)
h = uicontrol(Style = style, ...
String = string, ...
Units = "normalized", ...
Position = position, ...
Callback = callback, ...
Value = value, ...
Userdata = userdata, ...
Fontsize = 5 + 9*figsize, ...
FontWeight = "bold", ...
BackgroundColor = "w");
end
function [y,z,r2] = poly(h,t,p,x,w)
d = h.deg.UserData;
s = (t-1960)/60;
c = polyfit(s,p,d);
u = polyval(c,s);
s = (x-1960)/60;
y = polyval(c,s);
s = (w-1960)/60;
z = polyval(c,s);
r2 = 1-norm(p-u)^2/norm(p-mean(p))^2;
end
function [y,z] = hermite(model,t,p,x,w)
switch model
case 'pchip'
y = pchip(t,p,x);
z = pchip(t,p,w);
case 'spline'
y = spline(t,p,x);
z = spline(t,p,w);
case 'makima'
y = makima(t,p,x);
z = makima(t,p,w);
end
end
function [y,z,r2] = expo(model,t,p,x,w)
switch model
case 'exponential'
c = polyfit(log(t),log(p),1);
y = exp(polyval(c,log(x)));
z = exp(polyval(c,log(w)));
u = exp(polyval(c,log(t)));
case 'logistic'
logistic = @(a,b,c,t) a./(1+exp(-b*(t-c)));
fun = @(a) logistic_fit(a,t,p);
a = fminbnd(fun,500,1000);
[~,b,c] = logistic_fit(a,t,p);
y = logistic(a,b,c,x);
z = logistic(a,b,c,w);
u = logistic(a,b,c,t);
case 'gompertz'
gompertz = @(a,b,c,t) a*exp(-b*exp(-c*t));
fun = @(a)gompertz_fit(a,t,p);
a = fminbnd(fun,1000,6000);
[~,b,c] = gompertz_fit(a,t,p);
y = gompertz(a,b,c,x);
z = gompertz(a,b,c,w);
u = gompertz(a,b,c,t);
end
r2 = 1-norm(p-u)^2/norm(p-mean(p))^2;
end
function [rnorm,b,c] = logistic_fit(a,t,p)
% LOGISTIC_FIT. Objective function for one dimensional search.
% [rnorm,b,c] = logistic_fit(a,t,p)
% fits p with a./(1+exp(-b*(t-c))).
% Returns norm(fit - p).
% Uses linear fit of log(a./p-1) by -b*(t-c).
% Call a = fminbnd(@(a)logistic_fit(a,t,p),a_1,a_2).
logistic = @(a,b,c,t) a./(1+exp(-b*(t-c)));
A = [t ones(size(t))];
rhs = log(a./p-1);
f = A\rhs;
b = -f(1);
c = f(2)/b;
r = logistic(a,b,c,t) - p;
rnorm = norm(r);
end
function [rnorm,b,c] = gompertz_fit(a,t,p)
%GOMPERTZ_FIT. Objective function for one dimensional search.
% [rnorm,b,c] = gompertz_fit(a,t,p)
% fits p with a*exp(-b*exp(-c*t))).
% Returns norm(fit - p).
% Uses linear fit of log(log(a/p)) by log(-b)-c*t.
% Call a = fminbnd(@(a)gompertz_fit(a,t,p),a_1,a_2).
gompertz = @(a,b,c,t) a*exp(-b*exp(-c*t));
A = [t ones(size(t))];
rhs = log(log(a./p));
f = A\rhs;
b = exp(f(2));
c = -f(1);
r = gompertz(a,b,c,t) - p;
rnorm = norm(r);
end
function model_cb(~,~)
m = get(h.model,'value');
censusapp2024(models{m})
end
function predict_cb(arg,~)
sig = 2*(get(arg,'string') == '+') - 1;
w = get(h.predict,'userdata');
w = w + sig*(1 + 9*(w > 2039));
axis([1890 max(2030,w+20) 0 max(450,3*(w-1890))])
h.predict.String = sprintf('predict = %d',w);
h.predict.UserData = w;
model_cb
end
function noise_cb(arg,~)
sig = 2*(get(arg,'string') == '+') - 1;
level = get(h.noise,'userdata');
level = max(0,level + 10*sig);
h.noise.String = sprintf('noise = %d',level);
h.noise.UserData = level;
model_cb
end
function degree_cb(arg,~)
sig = 2*(get(arg,'string') == '+') - 1;
d = get(h.deg,'userdata');
d = min(12,max(0,d + sig));
h.deg.String = sprintf('degree = %d',d);
h.deg.UserData = d;
model_cb
end
function r2_cb(arg,~)
onoff = {'on','off'};
h.r2txt.Visible = onoff{2-get(arg,'value')};
end
function errs_cb(~,~)
model_cb
end
function info_cb(~,~)
web(['https://blogs.mathworks.com/cleve/2024/05/19/' ...
'a-sixty-year-old-program-for-predicting-the-future'], ...
'-notoolbar');
end
function help_cb(~,~)
helpwin('censusapp2024')
end
function restart_cb(~,~)
censusapp2024
end
function errest = errorestimates(model,t,p,x,y)
% Provide error estimates for censusapp2024
switch model
case 'polynomial'
d = get(h.deg,'userdata');
if d > 0
V(:,d+1) = ones(size(t));
s = (t-1960)/60;
for j = d:-1:1
V(:,j) = s.*V(:,j+1);
end
[Q,R] = qr(V);
r = p - V*(R\(Q'*p));
R = R(1:d+1,:);
RI = inv(R);
E = zeros(length(x),d+1);
s = (x-1960)/60;
for j = 1:d+1
E(:,j) = polyval(RI(:,j),s);
end
rho = norm(r,inf)+1;
e = rho*sqrt(1+diag(E*E'));
errest = [y-e; NaN; y+e];
else
errest = [y-NaN; NaN; y+NaN];
end
case {'exponential','logistic','gompertz'}
V = [ones(size(t)) log(t)];
[Q,R] = qr(V);
q = R\(Q'*log(p));
r = log(p) - V*q;
E = [ones(size(x)) log(x)]/R(1:2,1:2);
rho = norm(r,inf);
e = rho*sqrt(1+diag(E*E'));
errest = [y.*exp(-e); NaN; y.*exp(e)];
case {'pchip','spline','makima'}
n = length(t);
I = eye(n,n);
E = zeros(length(x),n);
for j = 1:n
switch model
case 'pchip'
E(:,j) = pchip(t,I(:,j),x);
case 'spline'
E(:,j) = spline(t,I(:,j),x);
case 'makima'
E(:,j) = makima(t,I(:,j),x);
end
end
rho = 1;
e = rho*sqrt(1+diag(E*E'));
errest = [y-e; NaN; y+e];
otherwise
errest = [y; NaN; y];
end % switch
end % errorestimates
end % censusapp2024