public
Last active

  • Download Gist
mortality.matlab
Matlab
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
%% 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);

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.