Systems-Modeling-Spring-2015-Notebooks/Thursday Apr 9 - Zombie SIR Model Fitting.py at master · bblais/Systems-Modeling-Spring-2015-Notebooks · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

from pyndamics import *
from pyndamics.emcee import *
import pandas

# <markdowncell>

# ## Data from homework

# <codecell>

with open('zombie_data.csv','w') as fid:
    fid.write("""time    zombies
    0.00    0
    3.00    1
    5.00    2
    6.00    2
    8.00    3
    10.00   3
    22.00   4
    22.20   6
    22.50   2
    24.00   3
    25.50   5
    26.00   12
    26.50   15
    27.50   25
    27.75   37
    28.50   25
    29.00   65
    29.50   80
    31.50   100""")

# <codecell>

data=pandas.read_csv('zombie_data.csv', delim_whitespace=True)
data

# <codecell>

sim=Simulation()
sim.add("S'=-beta*S*Z",300,plot=1)
sim.add("Z'=+beta*S*Z-gamma*Z",0.1,plot=2)
sim.add("R'=gamma*Z",0,plot=1)
sim.add_data(t=data['time'],Z=data['zombies'],plot=2)
sim.params(beta=.001,gamma=.03)
sim.run(0,50)

# <codecell>

model=MCMCModel(sim,beta=Uniform(0,0.001),gamma=Uniform(0,.1))

# <codecell>

model.run_mcmc(500)
model.plot_chains()

# <codecell>

model.set_initial_values('samples')
model.run_mcmc(500)
model.plot_chains()

# <codecell>

model.plot_distributions()

# <codecell>

model.triangle_plot()

# <codecell>

sim.run(0,50)

# <codecell>

sim.beta

# <codecell>

model.median_values

# <codecell>

model=MCMCModel(sim,beta=Uniform(0,0.01),gamma=Uniform(0,.1),initial_Z=Uniform(0,1))
model.run_mcmc(500)
model.plot_chains()

# <codecell>

model.set_initial_values('samples')
model.run_mcmc(500)
model.plot_chains()

# <codecell>

model.plot_distributions()

# <codecell>

model.triangle_plot()

# <codecell>

sim.noplots=True  # turn off the simulation plots
for i in range(500):
    model.draw()
    sim.run(0,50)
    plot(sim.t,sim.Z,'g-',alpha=.05)
sim.noplots=False  # gotta love a double-negative
plot(data['time'],data['zombies'],'bo')  # plot the data

# <codecell>