chap14soln


Modeling and Simulation in Python

Chapter 14

Copyright 2017 Allen Downey

License: Creative Commons Attribution 4.0 International


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 *

Code from previous chapters


In [2]:
def make_system(beta, gamma):
    """Make a system object for the SIR model.
    
    beta: contact rate in days
    gamma: recovery rate in days
    
    returns: System object
    """
    init = State(S=89, I=1, R=0)
    init /= np.sum(init)

    t0 = 0
    t_end = 7 * 14

    return System(init=init, t0=t0, t_end=t_end,
                  beta=beta, gamma=gamma)

In [3]:
def update_func(state, t, system):
    """Update the SIR model.
    
    state: State (s, i, r)
    t: time
    system: System object
    
    returns: State (sir)
    """
    s, i, r = state

    infected = system.beta * i * s    
    recovered = system.gamma * i
    
    s -= infected
    i += infected - recovered
    r += recovered
    
    return State(S=s, I=i, R=r)

In [4]:
def run_simulation(system, update_func):
    """Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    """
    init, t0, t_end = system.init, system.t0, system.t_end
    
    frame = TimeFrame(columns=init.index)
    frame.row[t0] = init
    
    for t in linrange(t0, t_end):
        frame.row[t+1] = update_func(frame.row[t], t, system)
    
    return frame

In [5]:
def calc_total_infected(results):
    """Fraction of population infected during the simulation.
    
    results: DataFrame with columns S, I, R
    
    returns: fraction of population
    """
    return get_first_value(results.S) - get_last_value(results.S)

In [6]:
def sweep_beta(beta_array, gamma):
    """Sweep a range of values for beta.
    
    beta_array: array of beta values
    gamma: recovery rate
    
    returns: SweepSeries that maps from beta to total infected
    """
    sweep = SweepSeries()
    for beta in beta_array:
        system = make_system(beta, gamma)
        results = run_simulation(system, update_func)
        sweep[system.beta] = calc_total_infected(results)
    return sweep

In [7]:
def sweep_parameters(beta_array, gamma_array):
    """Sweep a range of values for beta and gamma.
    
    beta_array: array of infection rates
    gamma_array: array of recovery rates
    
    returns: SweepFrame with one row for each beta
             and one column for each gamma
    """
    frame = SweepFrame(columns=gamma_array)
    for gamma in gamma_array:
        frame[gamma] = sweep_beta(beta_array, gamma)
    return frame

Contact number

Here's the SweepFrame from the previous chapter, with one row for each value of beta and one column for each value of gamma.


In [8]:
beta_array = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 , 1.1]
gamma_array = [0.2, 0.4, 0.6, 0.8]
frame = sweep_parameters(beta_array, gamma_array)
frame.head()


Out[8]:
0.2 0.4 0.6 0.8
0.1 0.010756 0.003642 0.002191 0.001567
0.2 0.118984 0.010763 0.005447 0.003644
0.3 0.589095 0.030185 0.010771 0.006526
0.4 0.801339 0.131563 0.020917 0.010780
0.5 0.896577 0.396409 0.046140 0.017640

In [9]:
frame.shape


Out[9]:
(11, 4)

The following loop shows how we can loop through the columns and rows of the SweepFrame. With 11 rows and 4 columns, there are 44 elements.


In [10]:
for gamma in frame.columns:
    column = frame[gamma]
    for beta in column.index:
        frac_infected = column[beta]
        print(beta, gamma, frac_infected)


0.1 0.2 0.010756340768063644
0.2 0.2 0.11898421353185373
0.3 0.2 0.5890954199973404
0.4 0.2 0.8013385277185551
0.5 0.2 0.8965769637207062
0.6 0.2 0.942929291399791
0.7 0.2 0.966299311298026
0.8 0.2 0.9781518959989762
0.9 0.2 0.9840568957948106
1.0 0.2 0.9868823507202488
1.1 0.2 0.988148177093735
0.1 0.4 0.0036416926514175607
0.2 0.4 0.010763463373360094
0.3 0.4 0.030184952469116566
0.4 0.4 0.131562924303259
0.5 0.4 0.3964094037932606
0.6 0.4 0.5979016626615987
0.7 0.4 0.7284704154876106
0.8 0.4 0.8144604459153759
0.9 0.4 0.8722697237137128
1.0 0.4 0.9116692168795855
1.1 0.4 0.9386802509510287
0.1 0.6 0.002190722188881611
0.2 0.6 0.005446688837466351
0.3 0.6 0.010771139974975585
0.4 0.6 0.020916599304195316
0.5 0.6 0.04614035896610047
0.6 0.6 0.13288938996079536
0.7 0.6 0.3118432512847451
0.8 0.6 0.47832565854255393
0.9 0.6 0.605687582114665
1.0 0.6 0.7014254793376209
1.1 0.6 0.7738176405451065
0.1 0.8 0.0015665254038139675
0.2 0.8 0.003643953969662994
0.3 0.8 0.006526163529085194
0.4 0.8 0.010779807499500693
0.5 0.8 0.017639902596349066
0.6 0.8 0.030291868201986594
0.7 0.8 0.05882382948158804
0.8 0.8 0.13358889291095588
0.9 0.8 0.2668895539427739
1.0 0.8 0.40375121210421994
1.1 0.8 0.519583469821867

Now we can wrap that loop in a function and plot the results. For each element of the SweepFrame, we have beta, gamma, and frac_infected, and we plot beta/gamma on the x-axis and frac_infected on the y-axis.


In [11]:
def plot_sweep_frame(frame):
    """Plot the values from a SweepFrame.
    
    For each (beta, gamma), compute the contact number,
    beta/gamma
    
    frame: SweepFrame with one row per beta, one column per gamma
    """
    for gamma in frame.columns:
        column = frame[gamma]
        for beta in column.index:
            frac_infected = column[beta]
            plot(beta/gamma, frac_infected, 'ro')

Here's what it looks like:


In [12]:
plot_sweep_frame(frame)

decorate(xlabel='Contact number (beta/gamma)',
         ylabel='Fraction infected')

savefig('figs/chap14-fig01.pdf')


Saving figure to file figs/chap14-fig01.pdf

It turns out that the ratio beta/gamma, called the "contact number" is sufficient to predict the total number of infections; we don't have to know beta and gamma separately.

We can see that in the previous plot: when we plot the fraction infected versus the contact number, the results fall close to a curve.

Analysis

In the book we figured out the relationship between $c$ and $s_{\infty}$ analytically. Now we can compute it for a range of values:


In [13]:
s_inf_array = linspace(0.0001, 0.9999, 101);

In [14]:
c_array = log(s_inf_array) / (s_inf_array - 1);

total_infected is the change in $s$ from the beginning to the end.


In [15]:
frac_infected = 1 - s_inf_array
frac_infected_series = Series(frac_infected, index=c_array);

Now we can plot the analytic results and compare them to the simulations.


In [16]:
plot_sweep_frame(frame)
plot(frac_infected_series, label='Analysis')

decorate(xlabel='Contact number (c)',
         ylabel='Fraction infected')

savefig('figs/chap14-fig02.pdf')


Saving figure to file figs/chap14-fig02.pdf

The agreement is generally good, except for values of c less than 1.

Exercises

Exercise: If we didn't know about contact numbers, we might have explored other possibilities, like the difference between beta and gamma, rather than their ratio.

Write a version of plot_sweep_frame, called plot_sweep_frame_difference, that plots the fraction infected versus the difference beta-gamma.

What do the results look like, and what does that imply?


In [17]:
# Solution

def plot_sweep_frame_difference(frame):
    for gamma in frame.columns:
        column = frame[gamma]
        for beta in column.index:
            frac_infected = column[beta]
            plot(beta - gamma, frac_infected, 'ro')

In [18]:
# Solution

plot_sweep_frame_difference(frame)

decorate(xlabel='Excess infection rate (infections-recoveries per day)',
         ylabel='Fraction infected',
         legend=False)



In [19]:
# Solution

# The results don't fall on a line, which means that if we know the difference between
# `beta` and `gamma`, but not their ratio, that's not enough to predict the fraction infected.

Exercise: Suppose you run a survey at the end of the semester and find that 26% of students had the Freshman Plague at some point.

What is your best estimate of c?

Hint: if you print frac_infected_series, you can read off the answer.


In [20]:
# Solution

frac_infected_series


Out[20]:
9.211261    0.999900
4.642296    0.989902
3.987365    0.979904
3.612133    0.969906
3.350924    0.959908
3.151808    0.949910
2.991711    0.939912
2.858363    0.929914
2.744467    0.919916
2.645332    0.909918
2.557767    0.899920
2.479505    0.889922
2.408879    0.879924
2.344627    0.869926
2.285771    0.859928
2.231541    0.849930
2.181315    0.839932
2.134590    0.829934
2.090947    0.819936
2.050040    0.809938
2.011573    0.799940
1.975299    0.789942
1.941002    0.779944
1.908499    0.769946
1.877628    0.759948
1.848249    0.749950
1.820238    0.739952
1.793487    0.729954
1.767898    0.719956
1.743384    0.709958
              ...   
1.181034    0.290042
1.173263    0.280044
1.165630    0.270046
1.158132    0.260048
1.150765    0.250050
1.143524    0.240052
1.136407    0.230054
1.129409    0.220056
1.122527    0.210058
1.115758    0.200060
1.109099    0.190062
1.102547    0.180064
1.096099    0.170066
1.089751    0.160068
1.083503    0.150070
1.077350    0.140072
1.071291    0.130074
1.065323    0.120076
1.059444    0.110078
1.053651    0.100080
1.047943    0.090082
1.042317    0.080084
1.036772    0.070086
1.031305    0.060088
1.025914    0.050090
1.020598    0.040092
1.015356    0.030094
1.010185    0.020096
1.005083    0.010098
1.000050    0.000100
Length: 101, dtype: float64

In [21]:
# Solution

# It looks like the fraction infected is 0.26 when the contact number is about 1.16

In [ ]: