In [66]:

    

%matplotlib inline

import scipy.integrate as spi
import numpy as np
import pylab as pl


    

In [81]:

    

#Person and Zombie

z_infecting = 0.4
p_killedby_z = 0.1
z_creates_i = 0.25
z_killedby_p = 0.75

#Person and Immune

i_infecting = 0.25
i_creates_i = 0.075

#Zombie no mates
z_starves = 0.05

#Immune and Zombie
i_killedby_z = 0.5
z_killedby_i = 0.2475


    

In [85]:

    

P = 0.8
Z = 0.0
I = 0.2
D = 0.0

if sum([P,Z,I,D]) != 1.0:
    print "Your proportions are wrong"

start_values = (P,Z,I,D)


    

In [69]:

    

t_start = 0.0
t_end = 100
t_step = 1


    

In [70]:

    

def time_range(t_start, t_end, t_step):
    '''
    Return step range of sampling points for iteration.
    '''
    t_ranges = np.arange(t_start, t_end + t_step, t_step)
    return t_ranges


    

In [71]:

    

t_ranges = time_range(t_start, t_end, t_step)


    

In [72]:

    

def calc_P(P,Z,I,D):
    out_P = - (z_infecting * (P * Z)) - (i_infecting * (P * I)) - (p_killedby_z * (P * Z)) - (z_creates_i * (P * Z)) - (i_creates_i * (P * I))
    return out_P


    

In [73]:

    

def calc_Z(P,Z,I,D):
    out_Z = (z_infecting * (P * Z)) + (i_infecting * (P * I)) - (z_killedby_p * (P * Z)) - (z_killedby_i * (I * Z)) - (z_starves * Z)
    return out_Z


    

In [74]:

    

def calc_I(P,Z,I,D):
    out_I = (z_creates_i * (P * Z)) + (i_creates_i * (P * I)) - (i_killedby_z * (I * Z))
    return out_I


    

In [75]:

    

def calc_D(P,Z,I,D):
    out_D = (p_killedby_z * (P * Z)) + (i_killedby_z * (I * Z)) + (z_killedby_p * (P * Z)) + (z_killedby_i * (Z * I)) + (z_starves * Z)
    return out_D


    

In [76]:

    

def diffeq(INP, t):
    OUT = np.zeros((4))
    mP,mZ,mI,mD = INP
    OUT[0] = calc_P(mP,mZ,mI,mD)
    OUT[1] = calc_Z(mP,mZ,mI,mD)
    OUT[2] = calc_I(mP,mZ,mI,mD)
    OUT[3] = calc_D(mP,mZ,mI,mD)
    return OUT


    

In [86]:

    

RES = spi.odeint(diffeq,start_values,t_ranges)


    

In [78]:

    

def plot_graph(RES):
    pl.close()
    pl.plot(RES[:,0], '-bo', label = 'People')
    pl.plot(RES[:,1], '-ro', label = 'Zombies')
    pl.plot(RES[:,2], '-go', label = 'Immune')
    pl.plot(RES[:,3], '-yo', label = 'Dead')
    pl.legend(loc = 0)
    pl.title('Zombie apocalypse')
    pl.xlabel('Time')
    pl.ylabel('Numbers')
    #pl.savefig('2.5_SIS-high.png', dpi = 900)
    pl.show()


    

In [79]:

    

observed = RES[:,0]
expected = [0.7,0.57582525,0.44145637,0.31872655,0.22099319,0.15023086,0.10178351,0.06947228,0.04806533,0.03381077]

total = 0

for i,item in enumerate(expected):
    total += (item - observed[i]) ** 2
    
print total


    

    




0.0332838936197

In [87]:

    

plot_graph(RES)plot_graph(RES)