Something along the line of:
s0, e0, i0 = 100., 50., 25.
st0, et0, it0 = [theano.shared(x) for x in [s0, e0, i0]]
C = np.array([3, 5, 8, 13, 21, 26, 10, 3], dtype=np.float64)
D = np.array([1, 2, 3, 7, 9, 11, 5, 1], dtype=np.float64)
def seir_one_step(st0, et0, it0, beta, gamma, delta):
bt0 = st0 * beta
ct0 = et0 * gamma
dt0 = it0 * delta
st1 = st0 - bt0
et1 = et0 + bt0 - ct0
it1 = it0 + ct0 - dt0
return st1, et1, it1
with pm.Model() as model:
beta = pm.Beta('beta', 2, 10)
gamma = pm.Beta('rho', 2, 10)
delta = pm.Beta('gamma', 2, 10)
(st, et, it), updates = theano.scan(
fn=seir_one_step,
outputs_info=[st0, et0, it0],
non_sequences=[beta, gamma, delta],
n_steps=len(C))
ct = pm.Binomial('c_t', et, gamma, observed=C)
dt = pm.Binomial('d_t', it, delta, observed=D)
trace = pm.sample()
with model:
bt = pm.Binomial('b_t', st, beta, shape=len(C))
ppc_trace = pm.sample_posterior_predictive(trace, var_names=['b_t'])
If you are interested in SEIR like model, Prediction of Danish Covid19 cases is a great example.