I have the following time series model

```
S[t+1] = S[t] - B[t]
E[t+1] = E[t] +B[t] - C[t]
I[t+1] = I[t+1] + C[t] - D[t]
B[t] ~ Binom(S[t], beta)
C[t] ~ Binom(E[t], gamma)
D[t] ~ Binom(I[t], delta)
beta, gamma, delta ~ Beta(a, b)
```

Only C and D are observed and we know S[0], E[0], I[0], a and b. What is the best way to represent this model in PyMC3?

Currently I do the following:

```
with pm.Model() as model:
s0, e0, i0 = 100, 50, 25
C = np.array([3, 5, 8, 13, 21, 26, 10, 3])
D = np.array([1, 2, 3, 7, 9, 11, 5, 1])
beta = pm.Beta('beta', 2, 10)
gamma = pm.Beta('rho', 2, 10)
delta = pm.Beta('gamma', 2, 10)
B = []
st, et, it = s0, e0, i0
for t in range(len(C)):
bt = pm.Binomial(f'b_{t}', st, beta)
ct = pm.Binomial(f'c_{t}', et, gamma, observed=C[t])
dt = pm.Binomial(f'd_{t}', it, delta, observed=D[t])
st = st - bt
et = et + bt - ct
it = it + ct - dt
B.append(bt)
step1 = pm.Metropolis(B)
step2 = pm.NUTS([beta, rho, gamma])
trace = pm.sample(1000, tune=1000, step=[step1, step2])
```

However, this creates one Metropolis sampler for each B[t] and does not scale well for a long time series. I would like to sample the entire B at once with one sampler. I also have a transition distribution for B in mind that might help the convergence.

How can I do it with PyMC3?