In [1]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *
from pandas import read_html
In [2]:
filename = 'data/World_population_estimates.html'
tables = read_html(filename, header=0, index_col=0, decimal='M')
table2 = tables[2]
table2.columns = ['census', 'prb', 'un', 'maddison',
'hyde', 'tanton', 'biraben', 'mj',
'thomlinson', 'durand', 'clark']
In [3]:
un = table2.un / 1e9
un.head()
In [4]:
census = table2.census / 1e9
census.head()
In [5]:
t_0 = get_first_label(census)
t_end = get_last_label(census)
elapsed_time = t_end - t_0
p_0 = get_first_value(census)
p_end = get_last_value(census)
total_growth = p_end - p_0
annual_growth = total_growth / elapsed_time
We can rewrite the code from the previous chapter using system objects.
In [6]:
system = System(t_0=t_0,
t_end=t_end,
p_0=p_0,
annual_growth=annual_growth)
And we can encapsulate the code that runs the model in a function.
In [7]:
def run_simulation1(system):
"""Runs the constant growth model.
system: System object
returns: TimeSeries
"""
results = TimeSeries()
results[system.t_0] = system.p_0
for t in linrange(system.t_0, system.t_end):
results[t+1] = results[t] + system.annual_growth
return results
We can also encapsulate the code that plots the results.
In [8]:
def plot_results(census, un, timeseries, title):
"""Plot the estimates and the model.
census: TimeSeries of population estimates
un: TimeSeries of population estimates
timeseries: TimeSeries of simulation results
title: string
"""
plot(census, ':', label='US Census')
plot(un, '--', label='UN DESA')
plot(timeseries, color='gray', label='model')
decorate(xlabel='Year',
ylabel='World population (billion)',
title=title)
Here's how we run it.
In [9]:
results = run_simulation1(system)
plot_results(census, un, results, 'Constant growth model')
Here's a more realistic model where the number of births and deaths is proportional to the current population.
In [10]:
def run_simulation2(system):
"""Run a model with proportional birth and death.
system: System object
returns: TimeSeries
"""
results = TimeSeries()
results[system.t_0] = system.p_0
for t in linrange(system.t_0, system.t_end):
births = system.birth_rate * results[t]
deaths = system.death_rate * results[t]
results[t+1] = results[t] + births - deaths
return results
I picked a death rate that seemed reasonable and then adjusted the birth rate to fit the data.
In [11]:
system.death_rate = 0.01
system.birth_rate = 0.027
Here's what it looks like.
In [12]:
results = run_simulation2(system)
plot_results(census, un, results, 'Proportional model')
savefig('figs/chap06-fig01.pdf')
The model fits the data pretty well for the first 20 years, but not so well after that.
run_simulation1 and run_simulation2 are nearly identical except the body of the loop. So we can factor that part out into a function.
In [13]:
def update_func1(pop, t, system):
"""Compute the population next year.
pop: current population
t: current year
system: system object containing parameters of the model
returns: population next year
"""
births = system.birth_rate * pop
deaths = system.death_rate * pop
return pop + births - deaths
The name update_func refers to a function object.
In [14]:
update_func1
Which we can confirm by checking its type.
In [15]:
type(update_func1)
run_simulation takes the update function as a parameter and calls it just like any other function.
In [16]:
def run_simulation(system, update_func):
"""Simulate the system using any update function.
system: System object
update_func: function that computes the population next year
returns: TimeSeries
"""
results = TimeSeries()
results[system.t_0] = system.p_0
for t in linrange(system.t_0, system.t_end):
results[t+1] = update_func(results[t], t, system)
return results
Here's how we use it.
In [17]:
t_0 = get_first_label(census)
t_end = get_last_label(census)
p_0 = census[t_0]
system = System(t_0=t_0,
t_end=t_end,
p_0=p_0,
birth_rate=0.027,
death_rate=0.01)
In [18]:
results = run_simulation(system, update_func1)
plot_results(census, un, results, 'Proportional model, factored')
Remember not to put parentheses after update_func1. What happens if you try?
Exercise: When you run run_simulation, it runs update_func1 once for each year between t_0 and t_end. To see that for yourself, add a print statement at the beginning of update_func1 that prints the values of t and pop, then run run_simulation again.
Since births and deaths get added up, we don't have to compute them separately. We can combine the birth and death rates into a single net growth rate.
In [19]:
def update_func2(pop, t, system):
"""Compute the population next year.
pop: current population
t: current year
system: system object containing parameters of the model
returns: population next year
"""
net_growth = system.alpha * pop
return pop + net_growth
Here's how it works:
In [20]:
system.alpha = system.birth_rate - system.death_rate
results = run_simulation(system, update_func2)
plot_results(census, un, results, 'Proportional model, combined birth and death')
Exercise: Maybe the reason the proportional model doesn't work very well is that the growth rate, alpha, is changing over time. So let's try a model with different growth rates before and after 1980 (as an arbitrary choice).
Write an update function that takes pop, t, and system as parameters. The system object, system, should contain two parameters: the growth rate before 1980, alpha1, and the growth rate after 1980, alpha2. It should use t to determine which growth rate to use. Note: Don't forget the return statement.
Test your function by calling it directly, then pass it to run_simulation. Plot the results. Adjust the parameters alpha1 and alpha2 to fit the data as well as you can.
In [21]:
# Solution goes here
In [22]:
# Solution goes here