Skip to content

Instantly share code, notes, and snippets.

@chris-taylor
Last active August 10, 2021 14:31
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save chris-taylor/5554669 to your computer and use it in GitHub Desktop.
%% Read in files and extrapolate mortality rate
x = csv.read('life.csv','format', '%f %f');
age = x{1};
prob = x{2};
model = stats.lm(log(prob(61:end)), log(age(61:100)));
pfun = @(a) util.if_(a<60, @()prob(a), util.if_(a>121, 1, exp(model.beta(1) + model.beta(2) * log(a))));
% model = stats.lm(log(prob(61:end)), age(61:end));
% pfun = @(a) util.if_(a<60, @()prob(a), util.if_(a>113, 1, exp(model.beta(1) + model.beta(2) * a)));
myprob = [];
for a = 1:122
myprob(a) = pfun(a);
end
%% Iterate the population distribution until convergence
P = 1e9; % total population of developed world
A = 122; % maximum age
T = 200; % number of years to run for
B = P/77.81; % birth rate
n = zeros(T, A); % population distribution in each generation
n(1,:) = repmat(B, 1, A); % uniform initial distribution
for t = 2:T
n(t,1) = B;
for a = 2:A
n(t,a) = (1-myprob(a-1)) * n(t-1,a-1);
end
end
%% Death simulation
pop = round(n(T, :)'); % Take the population to be the last generation from prev step
niters = 10000; % Number of years to simulate
nevents = 0; % Counter for number of times oldest person dies
ages = [];
clc
for t = 1:niters
new_pop = zeros(size(pop));
new_pop(1) = round(B); % New births
for a = 2:A
npeople = pop(a-1); % Number of people of age a-1 last year
pdeath = myprob(a-1); % Probability of death
if npeople * pdeath > 10 && npeople * (1-pdeath) > 10
% Normal approximation if appropriate
ndeaths = round(npeople*pdeath + sqrt(npeople*pdeath*(1-pdeath))*randn());
else
% Otherwise use binomial distribution
ndeaths = binornd(npeople, pdeath);
end
% Subtract the deaths from this year's cohort
new_pop(a) = pop(a-1) - ndeaths;
if new_pop(a) < 0 || isnan(new_pop(a))
keyboard
end
if ndeaths > 0 && all(pop(a:end) == 0) % is it possible that the oldest person has died?
nalive = pop(a-1);
for ii = 1:ndeaths
pOldestDied = 1 / nalive;
if rand() < pOldestDied
nevents = nevents + 1;
fprintf('Year is %d, oldest person has died age %d\n', t, a-1)
ages(end+1) = a-1;
end
nalive = nalive - 1;
end
end
end
pop = new_pop;
end
fprintf('There were %d events in %d years, for 1 every %.2f years\n', nevents, niters, niters/nevents);
[deathage,c] = util.count(ages);
bar(deathage,c);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment