pygom icon indicating copy to clipboard operation
pygom copied to clipboard

The passage of time is not a stochastic process: periodic functions may present some dificulty when using .simulate_jump()

Open m-d-grunnill opened this issue 3 years ago • 2 comments

I am trying to think of how you would stochastically simulate a model with a periodic term within pygom. The common_models module intialisis a determiistic SIS model with periodic transmision in the following way.

state = ['S', 'E', 'I', 'tau']
param_list = ['mu', 'alpha', 'gamma', 'beta_0', 'beta_1']
derived_param = [('beta_S', 'beta_0 * (1 + beta_1*cos(2*pi*tau))')]
ode = [
    Transition(origin='S', equation='mu - beta_S*S*I - mu*S',
               transition_type=TransitionType.ODE),
    Transition(origin='E', equation='beta_S*S*I - (mu + alpha)*E',
               transition_type=TransitionType.ODE),
    Transition(origin='I', equation='alpha*E - (mu + gamma)*I',
               transition_type=TransitionType.ODE),
    Transition(origin='tau', equation='1',
               transition_type=TransitionType.ODE)
    ]
# initialize the model
ode_obj = DeterministicOde(state, param_list,
                           derived_param=derived_param,
                           ode=ode)

From this example I would intialise a stochastic version of this model in the following way:

actual_states=['S', 'E', 'I', 'R']
state = actual_states + ['time']
param_list = ['mu', 'alpha', 'gamma', 'beta_0', 'beta_1']
derived_param = [('beta_S', 'beta_0 * (1 + beta_1*cos(2*pi*time))')]
transitions = [
    Transition(origin='S', destination='E', equation='beta_S*S*I',
               transition_type=TransitionType.T),
    Transition(origin='E', destination='I', equation='alpha*E',
               transition_type=TransitionType.T),
    Transition(origin='I', estination='R', equation='gamma*I',
               transition_type=TransitionType.T),
    ]
bd_list = [
   Transition(origin=states['S'][i], equation='mu*(S + E + I + R)', transition_type=TransitionType.B)]
   Transition(origin='time', equation='1', transition_type=TransitionType.B)
   ]

bd_list += [Transition(origin=item, equation='mu*'+item, transition_type=TransitionType.D) for item in actual_states]
    # initialize the model
    model = SimulateOde(state, param_list,
                               derived_param=derived_param,
                               transition=transition,
                               birth_death=bd_list)

The issue is that this would stochisticise time if the model.simulate_jump() method was called.

m-d-grunnill avatar Dec 07 '21 21:12 m-d-grunnill