<jupyter_start><jupyter_text>Modeling and Simulation in PythonChapter 10Copyright 2017 <NAME>License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>Under the hoodTo get a `DataFrame` and a `Series`, I'll read the world population data and select a column.`DataFrame` and `Series` contain a variable called `shape` that indicates the number of rows and columns.<jupyter_code>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']
table2.shape
census = table2.census / 1e9
census.shape
un = table2.un / 1e9
un.shape<jupyter_output><empty_output><jupyter_text>A `DataFrame` contains `index`, which labels the rows.  It is an `Int64Index`, which is similar to a NumPy array.<jupyter_code>table2.index<jupyter_output><empty_output><jupyter_text>And `columns`, which labels the columns.<jupyter_code>table2.columns<jupyter_output><empty_output><jupyter_text>And `values`, which is an array of values.<jupyter_code>table2.values<jupyter_output><empty_output><jupyter_text>A `Series` does not have `columns`, but it does have `name`.<jupyter_code>census.name<jupyter_output><empty_output><jupyter_text>It contains `values`, which is an array.<jupyter_code>census.values<jupyter_output><empty_output><jupyter_text>And it contains `index`:<jupyter_code>census.index<jupyter_output><empty_output><jupyter_text>If you ever wonder what kind of object a variable refers to, you can use the `type` function.  The result indicates what type the object is, and the module where that type is defined.`DataFrame`, `Int64Index`, `Index`, and `Series` are defined by Pandas.`ndarray` is defined by NumPy.<jupyter_code>type(table2)
type(table2.index)
type(table2.columns)
type(table2.values)
type(census)
type(census.index)
type(census.values)<jupyter_output><empty_output><jupyter_text>Optional exerciseThe following exercise provides a chance to practice what you have learned so far, and maybe develop a different growth model.  If you feel comfortable with what we have done so far, you might want to give it a try.**Optional Exercise:** On the Wikipedia page about world population estimates, the first table contains estimates for prehistoric populations.  The following cells process this table and plot some of the results. Select `tables[1]`, which is the second table on the page.<jupyter_code>table1 = tables[1]
table1.head()<jupyter_output><empty_output><jupyter_text>Not all agencies and researchers provided estimates for the same dates.  Again `NaN` is the special value that indicates missing data.<jupyter_code>table1.tail()<jupyter_output><empty_output><jupyter_text>Some of the estimates are in a form we can't read as numbers.  We could clean them up by hand, but for simplicity I'll replace any value that has an `M` in it with `NaN`.<jupyter_code>table1.replace('M', np.nan, regex=True, inplace=True)<jupyter_output><empty_output><jupyter_text>Again, we'll replace the long column names with more convenient abbreviations.<jupyter_code>table1.columns = ['prb', 'un', 'maddison', 'hyde', 'tanton', 
                  'biraben', 'mj', 'thomlinson', 'durand', 'clark']<jupyter_output><empty_output><jupyter_text>This function plots selected estimates.<jupyter_code>def plot_prehistory(table):
    """Plots population estimates.
    
    table: DataFrame
    """
    plot(table.prb, 'ro', label='PRB')
    plot(table.un, 'co', label='UN')
    plot(table.hyde, 'yo', label='HYDE')
    plot(table.tanton, 'go', label='Tanton')
    plot(table.biraben, 'bo', label='Biraben')
    plot(table.mj, 'mo', label='McEvedy & Jones')<jupyter_output><empty_output><jupyter_text>Here are the results.  Notice that we are working in millions now, not billions.<jupyter_code>plot_prehistory(table1)
decorate(xlabel='Year', 
         ylabel='World population (millions)',
         title='Prehistoric population estimates')<jupyter_output><empty_output><jupyter_text>We can use `xlim` to zoom in on everything after Year 0.<jupyter_code>plot_prehistory(table1)
decorate(xlim=[0, 2000], xlabel='Year', 
         ylabel='World population (millions)',
         title='Prehistoric population estimates')<jupyter_output><empty_output><jupyter_text>See if you can find a model that fits these data well from Year -1000 to 1940, or from Year 1 to 1940.How well does your best model predict actual population growth from 1950 to the present?<jupyter_code># Solution

def update_func_prop(pop, t, system):
    """Compute the population next year with proportional growth.
    
    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
# Solution

t_0 = 1
p_0 = table1.biraben[t_0]

prehistory = System(t_0=t_0, 
                    t_end=2016,
                    p_0=p_0,
                    alpha=0.0011)
# Solution

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
# Solution

results = run_simulation(prehistory, update_func_prop)
plot_prehistory(table1)
plot(results, color='gray', label='model')
decorate(xlim=[0, 2000], xlabel='Year', 
         ylabel='World population (millions)',
         title='Prehistoric population estimates')
# Solution

plot(census, ':', label='US Census')
plot(un, '--', label='UN DESA')
plot(results / 1000, color='gray', label='model')
decorate(xlim=[1950, 2016], xlabel='Year', 
         ylabel='World population (billions)',
         title='Prehistoric population estimates')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10721966205837173, 'fraction_numerical': 0.016129032258064516, 'mean_word_length': 5.062383612662942, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '13974', 'license': 'MIT', 'n_tokens_mistral': 2118, 'n_tokens_neox': 2019, 'n_words': 596, 'path': 'code/soln/chap10soln.ipynb', 'repo_name': 'arunkhattri/ModSimPy', 'size': '6516'}