# Modeling and Simulation in Python

Chapter 14



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

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)




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 [ ]: