I built a dual SIR model where my infectivity parameter varied sinusoidally over time so that I could simulate years of data rather than just one period. But I’m having trouble with the code.
Here is the code:
def dual_SIR_sunode(self, t, y, p):
# define β1/2(t) = B * (1 + b * sin(2 * pi * t / 365))
beta1_t = p.B1 * (1 + p.b1 * pm.math.sin(2 * np.pi * t / 52))
beta2_t = p.B2 * (1 + p.b2 * pm.math.sin(2 * np.pi * t / 52))
# I1_prime and I2_prime
I1_prime = y.IS + y.II + y.IR
I2_prime = y.SI + y.II + y.RI
return {
'SS': p.sigma1 * y.RS + p.sigma2 * y.SR - beta1_t * y.SS * I1_prime / y.n - beta2_t * y.SS * I2_prime / y.n,
'IS': beta1_t * y.SS * I1_prime / y.n + p.sigma2 * y.IR - p.gamma1 * y.IS - beta2_t * (1 - p.phi12) * y.IS * I2_prime / y.n,
'RS': p.sigma2 * y.RR + p.gamma1 * y.IS - p.sigma1 * y.RS - beta2_t * y.SS * I2_prime / y.n,
'SI': p.sigma1 * y.RI + beta2_t * y.SS * I2_prime / y.n - beta1_t * (1 - p.phi21) * y.SI * I1_prime / y.n - p.gamma2 * y.SI ,
'II': beta1_t * (1 - p.phi21) * y.SI * I1_prime / y.n + beta2_t * (1 - p.phi12) * y.IS * I2_prime / y.n - p.gamma1 * y.II - p.gamma2 * y.II,
'RI': p.gamma1 * y.II + beta2_t * y.RS * I2_prime / y.n - p.sigma1 * y.RI - p.gamma2 * y.RI,
'SR': p.sigma1 * y.RR + p.gamma2 * y.SI - beta1_t * y.SR * I1_prime / y.n - p.sigma2 * y.SR ,
'IR': beta1_t * y.SR * I1_prime / y.n + p.gamma2 * y.II - p.gamma1 * y.IR - p.sigma2 * y.IR ,
'RR': p.gamma1 * y.IR + p.gamma2 * y.RI - p.sigma1 * y.RR - p.sigma2 * y.RR
}
def run_SEIRS_model(self):
with pm.Model() as self.model:
# Priors for unknown model parameters
B1 = pm.Normal('B1', mu=0.5, sigma=0.1)
b1 = pm.Normal('b1', mu=0.1, sigma=0.02)
B2 = pm.Normal('B2', mu=0.5, sigma=0.1)
b2 = pm.Normal('b2', mu=0.1, sigma=0.02)
phi12 = pm.Beta('phi12', alpha=2, beta=5)
phi21 = pm.Beta('phi21', alpha=2, beta=5)
sigma1 = pm.Normal('sigma1', mu=0.1, sigma=0.05)
sigma2 = pm.Normal('sigma2', mu=0.1, sigma=0.05)
gamma1 = pm.Normal('gamma1', mu=0.1, sigma=0.02)
gamma2 = pm.Normal('gamma2', mu=0.1, sigma=0.02)
res, _, problem, solver, _, _ = sunode.wrappers.as_pytensor.solve_ivp(
y0={
# The initial conditions of the ode. Each variable
# needs to specify a theano or numpy variable and a shape.
# This dict can be nested.
'SS': (self.init_y['SS'], ()),
'IS': (self.init_y['IS'], ()),
'RS': (self.init_y['RS'], ()),
'SI': (self.init_y['SI'], ()),
'II': (self.init_y['II'], ()),
'RI': (self.init_y['RI'], ()),
'SR': (self.init_y['SR'], ()),
'IR': (self.init_y['IR'], ()),
'RR': (self.init_y['RR'], ()),
'n': (self.init_y['n'], ())
},
params={
# Each parameter of the ode. sunode will only compute derivatives
# with respect to theano variables. The shape needs to be specified
# as well. It it infered automatically for numpy variables.
# This dict can be nested.
'B1': (B1, ()),
'b1': (b1, ()),
'B2': (B2, ()),
'b2': (b2, ()),
'phi12': (phi12, ()),
'phi21': (phi21, ()),
'sigma1': (sigma1, ()),
'sigma2': (sigma2, ()),
'gamma1': (gamma1, ()),
'gamma2': (gamma2, ()),
'_dummy': (np.array(1.), ()),
},
# A functions that computes the right-hand-side of the ode using
# sympy variables.
rhs=self.dual_SIR_sunode,
# The time points where we want to access the solution
tvals=self.times[1:],
t0=self.times[0]
)
# The solution of the ode. This is a dictionary with the same structure
# as the initial conditions.
I1_total = res['IS'] + res['II'] + res['IR']
I2_total = res['SI'] + res['II'] + res['RI']
# Define the likelihood function
I1_obs = pm.Normal('I1_obs', mu=I1_total, sigma=0.1, observed=self.Dieases1_num[1:])
I2_obs = pm.Normal('I2_obs', mu=I2_total, sigma=0.1, observed=self.Dieases2_num[1:])
# Sample the model
with self.model:
trace = pm.sample(**self.samples_params)
return trace
However error happens on # define β1/2(t) = B * (1 + b * sin(2 * pi * t / 365)): NotImplementedError: Cannot convert 0.120830486676531time to a tensor variable.
I know this won’t work, but I can’t think of any other way to pass in the t parameter as a numerical value. I tried defind it outside the dual_SIR_sunode function then I can’t get the time parameter, like this:
Maybe someone has experience in this regard? Thanks to all who provided suggestions!
