Stochastic simulation

Examples taken from https://arxiv.org/pdf/1803.06934.pdf (see page 11 for stochastic simulations).

Examples are performed on an SIR model.

$\frac{dS}{dt} = -\beta S I $

$\frac{dI}{dt} = \beta S I - \gamma I$

$\frac{dR}{dt} = \gamma I$



In [1]:

    
import pygom
import pkg_resources
print('PyGOM version %s' %pkg_resources.get_distribution('pygom').version)









    



PyGOM version 0.1.3

Parameter stochasticity

Parameter values are sampled from a distribution. Deterministic solutions are found for a given set of parameter values.

In this example, $\beta$ and $\gamma$ are sampled from a gamma distribution.



In [2]:

    
from pygom import Transition, TransitionType, SimulateOde
import numpy as np

# construct model 
states = ['S', 'I', 'R']
params = ['beta', 'gamma', 'N']
transitions = [Transition(origin='S', destination='I', equation='beta*S*I/N', 
                          transition_type=TransitionType.T),
               Transition(origin='I', destination='R', equation='gamma*I', 
                          transition_type=TransitionType.T)]

model_p = SimulateOde(states, params, transition=transitions)



In [14]:

    
# initial conditions 
N = 7781984.0
in_inf = round(0.0000001*N)
init_state = [N - in_inf, in_inf, 0.0]

# time 
t = np.linspace (0 , 50 , 101)

# deterministic parameter values
param_evals = [('beta', 3.6), ('gamma', 0.2), ('N', N)]



In [4]:

    
# define parameter distributions
from pygom.utilR import rgamma

d = dict()
d['beta'] = (rgamma, {'shape':3600.0, 'rate':1000.0})
d['gamma'] = (rgamma, {'shape':1000.0, 'rate':500.0})
d['N'] = N



In [5]:

    
model_p.parameters = d
model_p.initial_values = (init_state, t[0])



In [6]:

    
# solve for 10 parameter sets 
Ymean, Yall = model_p.simulate_param(t[1::], iteration=10, full_output=True)



In [9]:

    
# plot solutions
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [12, 6]

fig, (ax1, ax2, ax3) = plt.subplots(1,3)
for i in range(np.shape(Yall)[0]):
    ax1.plot(t, Yall[i][:,0])
    ax2.plot(t, Yall[i][:,1])
    ax3.plot(t, Yall[i][:,2])
ax1.plot(t, Ymean[:,0], linewidth=3,color='k')
ax2.plot(t, Ymean[:,1], linewidth=3,color='k')
ax3.plot(t, Ymean[:,2], linewidth=3,color='k')
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax3.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax1.set_title('Susceptible')
ax2.set_title('Infected')
ax3.set_title('Removed')









    Out[9]:





Text(0.5, 1.0, 'Removed')



In [8]:

    
#!TODO Add legend with beta and gamma parameter

Jump process

Movements between states are discrete, where the probability of transition is given by

$Pr(process\ j\ jump\ within\ time\ \tau) = \lambda_j \exp^{-\lambda_j \tau}$.



In [16]:

    
# construct model
model_j = SimulateOde(states, params, transition=transitions)
model_j.parameters = param_evals
model_j.initial_values = (init_state, t[0])



In [ ]:

    
# run 10 simulations
simX, simT = model_j.simulate_jump(t[1::], iteration=10, full_output=True)









    



Parallel simulation
Revert to serial



In [ ]: