%  This routine computes a logistic meta-model for the U.S. census problem:
%  ------------------------------------------------------------------------
%
diary census.dry
echo on
%
%  Define the data
%
t = 1850:10:2000;
t = t';
p = [23.19,31.44,38.56,50.16,62.95,75.96,91.97,105.7,122.8,131.7,150.7,179.3,203.2,225.6,267.0,285.4];
p = 1E6*p';
%
%  Compute a decent estimate for the derivatives
%
[n,m] = size(p);
pdot = zeros(size(p));
for i=2:n-1,
  t1 = t(i-1:i+1);
  p1 = p(i-1:i+1);
  vdm = [ ones(3,1) , t1 , t1.*t1 ];
  coef = vdm\p1;
  pdot(i) = coef(2) + 2*coef(3)*t(i);
  if i==2,
    pdot(1) = coef(2) + 2*coef(3)*t(1);
  end,
  if i==n-1,
    pdot(n) = coef(2) + 2*coef(3)*t(n);
  end,
end
%
%  Compute the meta-model
%
vdm = [ p , p.*p ];
coef = vdm\pdot;
a = coef(1)
b = coef(2)
%
%  Compute the maximum value of the U.S. population
%
pmax = -coef(1)/coef(2)
%
%  Compute pm(t) using the meta-model
%
one = ones(size(t));
t0 = t(1);
pm0 = p(1);
c2 = 1/pmax;
c1 = 1/pm0 - c2;
pm = one ./ ( c1*exp(-coef(1)*(t - t0*one)) + c2*one );
%
%  Validate the meta-model by comparing true and fitted populations
%
plot(t,[p,pm])
grid on
title('Comparison of true and modelled US populations')
xlabel('time')
ylabel('population')
print -djpeg census.jpg
%
%  Find the year when the second derivative is zero
%
bpm0 = coef(2)*pm0;
tcrit = round(-log(-bpm0/(bpm0 + coef(1)))/coef(1)) + t0
%
diary off
return