Stochastic SIR-model with binomial likelihood

I’m currently working on this stochastic SIR-model:

\begin{align} \begin{split} \Delta S &= - \text{Po}(S_k \cdot p_I(t))\\ \Delta I &= \text{Po}(S_k \cdot p_I(t)) - \text{Bin}(I_k, p_R)\\ \Delta R &= \text{Bin}(I_k, p_R) \end{split} \end{align}
(Where the poisson distribution might aswell be swapped out with a binomial)
I generate data with a deterministic model:
\begin{align} \begin{split} \dot{S} &= S - \beta \frac{SI}{N}\\ \dot{I} &= I + \beta \frac{SI}{N} - \alpha I\\ \dot{R} &= R + \alpha I \end{split} \end{align}
Resulting in measurements y(t)=I(t). My goal here is to illustrate parameter uncertainty, and I’m planning to do it with the stochastic model. This results in two problems:

  1. DifferentialEquation(.) only works with continous models. How to I implement a custom, discrete differential equation?

There is a guide on ‘black box’ likelihood functions, which may solve this issue: Using a “black box” likelihood function — PyMC3 3.10.0 documentation

  1. How would I choose and define a likelihood function for my measurement in this case? The probability distributions are state-dependent.