chris-taylor · August 10, 2021 14:31 · May 10, 2013 · May 10, 2013 · May 10, 2013 · May 10, 2013
diff --git a/mortality.matlab b/mortality.matlab
@@ -18,7 +18,7 @@ end
 
 %% Iterate the population distribution until convergence
 
-P = 300e5;                  % total population of developed world
+P = 1e9;                  % total population of developed world
 
 A = 122;                  % maximum age
 T = 200;                  % number of years to run for

diff --git a/mortality.matlab b/mortality.matlab
@@ -7,17 +7,22 @@ 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 = 7e9;                  % this is the total population
+P = 300e5;                  % total population of developed world
 
 A = 122;                  % maximum age
 T = 200;                  % number of years to run for
-B = 1.286 * P;            % birth rate
+B = P/77.81;              % birth rate
 
 n = zeros(T, A);          % population distribution in each generation
 n(1,:) = repmat(B, 1, A); % uniform initial distribution
@@ -34,12 +39,13 @@ end
 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) = B; % New births
+    new_pop(1) = round(B); % New births
 
     for a = 2:A
         npeople = pop(a-1);    % Number of people of age a-1 last year
@@ -53,12 +59,17 @@ for t = 1:niters
         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
@@ -71,4 +82,5 @@ 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);
diff --git a/mortality.matlab b/mortality.matlab
@@ -54,15 +54,13 @@ for t = 1:niters
         % Subtract the deaths from this year's cohort
         new_pop(a) = pop(a-1) - ndeaths;
         if ndeaths > 0 && all(pop(a:end) == 0) % is it possible that the oldest person has died?
-            K = ndeaths;
             nalive = pop(a-1);
-            while K > 0
+            for ii = 1:ndeaths
                 pOldestDied = 1 / nalive;
                 if rand() < pOldestDied
                     nevents = nevents + 1;
                 end
-                K = K - 1;
-                nalive = nalive-1;
+                nalive = nalive - 1;
             end
         end
     end

diff --git a/mortality.matlab b/mortality.matlab
@@ -53,11 +53,16 @@ for t = 1:niters
         end
         % Subtract the deaths from this year's cohort
         new_pop(a) = pop(a-1) - ndeaths;
-        if ndeaths > 1 && all(pop(a:end) == 0) % is it possible that the oldest person has died?
-            pOldestDied = ndeaths / pop(a-1);  % probability that the oldest person has died
-            if rand() < pOldestDied
-                fprintf('Year is %d, oldest person in the world just died aged %d\n', t, a-1);
-                nevents = nevents + 1;
+        if ndeaths > 0 && all(pop(a:end) == 0) % is it possible that the oldest person has died?
+            K = ndeaths;
+            nalive = pop(a-1);
+            while K > 0
+                pOldestDied = 1 / nalive;
+                if rand() < pOldestDied
+                    nevents = nevents + 1;
+                end
+                K = K - 1;
+                nalive = nalive-1;
             end
         end
     end

diff --git a/mortality.matlab b/mortality.matlab
@@ -0,0 +1,71 @@
+%% 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))));
+for a = 1:122
+    myprob(a) = pfun(a);
+end
+
+%% Iterate the population distribution until convergence
+
+P = 7e9;                  % this is the total population
+
+A = 122;                  % maximum age
+T = 200;                  % number of years to run for
+B = 1.286 * P;            % 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
+clc
+
+for t = 1:niters
+
+    new_pop = zeros(size(pop));
+    new_pop(1) = 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 ndeaths > 1 && all(pop(a:end) == 0) % is it possible that the oldest person has died?
+            pOldestDied = ndeaths / pop(a-1);  % probability that the oldest person has died
+            if rand() < pOldestDied
+                fprintf('Year is %d, oldest person in the world just died aged %d\n', t, a-1);
+                nevents = nevents + 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);
+
+
No results found