Last active
August 10, 2021 14:31
-
-
Save chris-taylor/5554669 to your computer and use it in GitHub Desktop.
Code to answer this math.stackexchange question: http://math.stackexchange.com/questions/349155/how-often-does-it-happen-that-the-oldest-person-alive-dies/387581
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
%% 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