Instantly share code, notes, and snippets.

# chris-taylor/mortality.matlab Last active Dec 17, 2015

 %% 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);