Pyndamics provides a way to describe a dynamical system in terms of the differential equations, or the stock-flow formalism. It is a wrapper around the Scipy odeint function, with further functionality for time plots, phase plots, and vector fields. The MCMC component of this package uses pymc (version 2 right now): http://pymc-devs.github.io/pymc/. Page for this package: https://code.google.com/p/pyndamics/
from pyndamics import Simulation
from pyndamics.mcmc import MCMCModel
sim=Simulation()
sim.add("S'=-beta*S*I",1,plot=1)
sim.add("I'=beta*S*I-zeta*I",.001,plot=1)
sim.add("R'=zeta*I",0,plot=1)
sim.params(beta=5,zeta=1)
sim.run(0,10)
<matplotlib.figure.Figure at 0xccfe7f0>
sim=Simulation()
sim.add("S'=-beta*S*I",1,plot=1)
sim.add("E'=beta*S*I-zeta*E",0,plot=1)
sim.add("I'=zeta*E-alpha*I",.001,plot=1)
sim.add("R'=alpha*I",0,plot=1)
sim.params(alpha=.3,beta=10,zeta=.5)
sim.run(0,10)
<matplotlib.figure.Figure at 0xce19370>
Notice that no matter what the parameters are changed to, Z (zombies) always win.
sim=Simulation()
sim.add("S'=Pi-beta*S*Z-delta*S",500,plot=1) #S (Susceptible)
sim.add("Z'=beta*S*Z+zeta*R-alpha*S*Z",.002,plot=1) #Z (Zombie)
sim.add("R'=delta*S+alpha*S*Z-zeta*R",1,plot=False) #R (Removed)
sim.params(alpha=.005,beta=.0095,zeta=.05, delta=.01,Pi=0) #parameters changed to match the Munz et al. (2009) figures
sim.run(0,10)
<matplotlib.figure.Figure at 0xcb635f0>
Movie "data" from Night of the Living Dead
t=array([0,1,1.5,3,4.5,5,5.75,5.9,10])
zombies=array([1,1,3,8,10,20,28,30,40])
Simulation using data (we've plugged in the resulting estimated parameters)
sim=Simulation()
sim.add("S'=-beta*S*Z-delta*S",178.5,plot=1)
sim.add("E'=beta*S*Z-zeta*E",0,plot=False)
sim.add("Z'=zeta*E-alpha*S*Z",1,plot=1)
sim.add("R'=alpha*S*Z+delta*S",0,plot=False)
sim.params(alpha=.0342,beta=.0445,zeta=4.63, delta=0.0001)
sim.add_data(t=t,Z=zombies,plot=1)
sim.run(0,10)
<matplotlib.figure.Figure at 0xcdcfa70>
MCMC parameter estimation for the initial susceptible (S) population
model=MCMCModel(sim,{'initial_S':[50,500]})
model.fit(iter=100000)
model.plot_distributions()
[****************100%******************] 100000 of 100000 complete Time taken: 17 m, 59.59 s
MCMC parameter estimation for alpha (rate of zombies being permanently removed), beta (rate of susceptibles becoming infected), zeta (the rate of infected into becoming zombies), and delta (suicide rate among susceptibles)
#check priors
model=MCMCModel(sim,{'alpha':[0,5],'beta':[0,20],'zeta':[0,10]})
model.fit(iter=100000)
model.plot_distributions()
fig=sim.figures[0]
[****************100%******************] 100000 of 100000 complete Time taken: 1 h, 53 m, 40.28 s
plot(model['S'].value,'o')
plot(model['S'].stats()['mean'],'--')
plot(model['S'].stats()['95% HPD interval'],'r:')
plot(model['Z'].value,'o')
plot(model['Z'].stats()['mean'],'--')
plot(model['Z'].stats()['95% HPD interval'],'r:')
plot(model['R'].value,'o')
plot(model['R'].stats()['mean'],'--')
plot(model['R'].stats()['95% HPD interval'],'r:')
xlabel('Time')
ylabel('Population')
<matplotlib.text.Text at 0xd0dfcf0>
Joint distributions for alpha and Beta parameters - notice the direct correlation between the two terms
model.plot_joint_distribution('alpha','beta',show_prior=False)
Joint distributions for alpha and zeta parameters - notice that there is no co-dependency of the two terms
model.plot_joint_distribution('alpha','zeta',show_prior=False)
Movie "data" from Shaun of the Dead
t=array([0,3,5,6,8,10,22,22.2,22.5,24,25.5,26,26.5,27.5,27.75,28.5,29,29.5,31.5])
zombies=array([0,1,2,2,3,3,4,6,2,3,5,12,15,25,37,25,65,80,100])
Simulation using data (we've plugged in the resulting estimated parameters)
sim=Simulation()
sim.add("S'=-beta*S*Z",508.2,plot=1)
sim.add("E'=beta*S*Z-zeta*E",0,plot=0)
sim.add("Z'=zeta*E-alpha*S*Z",.000347759,plot=1)
sim.add("R'=alpha*S*Z",0,plot=False)
sim.params(alpha=2.96e-8,beta=0.000808995,zeta=60)
sim.add_data(t=t,Z=zombies,plot=1)
sim.run(0,50)
<matplotlib.figure.Figure at 0xd0f1b10>
MCMC parameter estimation for the initial susceptible (S) population
model=MCMCModel(sim,{'initial_S':[0,1000]})
model.fit(iter=100000)
model.plot_distributions()
[****************100%******************] 100000 of 100000 complete Time taken: 11 m, 27.98 s
MCMC parameter estimation for alpha (rate of zombies being permanently removed), beta (rate of susceptibles becoming infected), and zeta (the rate of infected into becoming zombies)
model=MCMCModel(sim,{'alpha':[0,5],'beta':[0,5], 'zeta':[0,10]})
model.fit(iter=100000)
model.plot_distributions()
[****************100%******************] 100000 of 100000 complete Time taken: 4 h, 14 m, 54.04 s
Data from GoogleTrends based on influenza records in Argentina (2003)
t=array([1,8,15,22,29,36,43,50,57,64,71,78,85,92,99,106,113,120,127,134,141,148,155,162,169,176,183,190,197,204,211,218,225,232,239,246,253,260,267,274,281,288,295,302,309,316,323,330,337,344])
flu=array([136,145,141,135,134,136,134,153,158,187,193,186,209,239,241,232,251,253,261,258,278,348,325,307,239,229,235,225,201,192,224,210,204,213,217,205,174,178,178,184,166,149,133,123,121,109,103,100,98,98])
Simulation using data (we've plugged in the resulting estimated parameters)
S0=10000
sim=Simulation()
sim.add("S'=-beta/S0*S*I",S0,plot=False)
sim.add("I'=beta/S0*S*I-alpha*I",136,plot=1)
sim.add("R'=alpha*I",0,plot=False)
sim.params(S0=S0, alpha=0.04334,beta=.05041)
sim.add_data(t=t,I=flu,plot=1)
sim.run(0,350)
<matplotlib.figure.Figure at 0xc96a090>
MCMC parameter estimation for alpha (rate of recovery) and beta (rate of infection)
model=MCMCModel(sim,{'alpha':[0,.3],'beta':[0,.3]})
model.fit(iter=100000)
model.plot_distributions()
[****************100%******************] 100000 of 100000 complete Time taken: 16 m, 49.99 s
Change Beta to be slightly smaller, where alpha > Beta (0.0393>.03)
sim=Simulation()
sim.add("S'=-beta*S*Z",180,plot=1)
sim.add("Z'=zeta*R-alpha*S*Z",1,plot=1)
sim.add("R'=beta*S*Z-zeta*R",0,plot=False)
sim.add("X'=alpha*S*Z",0,plot=False)
sim.params(alpha=.0393,beta=.03,zeta=6.15)
sim.run(0,10)
<matplotlib.figure.Figure at 0xd454f10>
Simulating military intervention in the film
def alpha(t):
if t>31.5: #t=31.5 because the zombie population reached about 100 or so at the end
return 0.0115
else:
return 0.00115
sim=Simulation()
sim.add("S'=-beta*S*Z",500,plot=1)
sim.add("Z'=beta*S*Z-alpha(t)*S*Z",.000347759,plot=1)
sim.add("X'=alpha(t)*S*Z",0,plot=False)
sim.params(beta=.00199)
sim.functions(alpha)
sim.run(0,50)
<matplotlib.figure.Figure at 0xd43f830>
Simulating military intervention had they come at 40 hours (8.5 hours later than in the movie)
def alpha(t):
if t>42: #t=40
return 0.0115
else:
return 0.00115
sim=Simulation()
sim.add("S'=-beta*S*Z",500,plot=1)
sim.add("Z'=beta*S*Z-alpha(t)*S*Z",.000347759,plot=1)
sim.add("X'=alpha(t)*S*Z",0,plot=False)
sim.params(beta=.00199)
sim.functions(alpha)
sim.run(0,50)
<matplotlib.figure.Figure at 0xd2f2ad0>