Estoy tratando de armar un modelo de un sistema dinámico en PyMC3, para inferir dos parámetros. El modelo es el SIR básico, comúnmente utilizado en epidemiología:

dS / dt = - r0 * g * S * I

dI / dt = g * I (r * S - 1)

donde r0 y g son parámetros a inferir. Hasta ahora, no puedo llegar muy lejos. Los únicos ejemplos que he visto de armar una cadena de Markov como ésta arrojan errores acerca de que la recursión es demasiado profunda. Aquí está mi código de ejemplo.

# Time
t = np.linspace(0, 8, 200)

# Simulated observation
def SIR(y, t, r0, gamma) :
    S = - r0 * gamma * y[0] * y[1]
    I = r0 * gamma * y[0] * y[1] - gamma * y[1]
    return [S, I]

# Currently no noise, we just want to infer params r0 = 16 and g = 0.5
solution = odeint(SIR, [0.99, 0.01, 0], t, args=(16., 0.5))


with pymc.Model() as model :
    r0 = pymc.Normal("r0", 15, sd=10)
    gamma = pymc.Uniform("gamma", 0.3, 1.)

    # Use forward Euler to solve
    dt = t[1] - t[0]

    # Initial conditions
    S = [0.99]
    I = [0.01]

    for i in range(1, len(t)) :
        S.append(pymc.Normal("S%i" % i, \
                         mu = S[-1] + dt * (-r0 * gamma * S[-1] * I[-1]), \
                         sd = solution[:, 0].std()))
        I.append(pymc.Normal("I%i" % i, \
                         mu = I[-1] + dt * ( r0 * gamma * S[-1] * I[-1] - gamma * I[-1]), \
                         sd = solution[:, 1].std()))

    Imcmc = pymc.Normal("Imcmc", mu = I, sd = solution[:, 1].std(), observed = solution[:, 1])

    #start = pymc.find_MAP()
    trace = pymc.sample(2000, pymc.NUTS())

Cualquier ayuda sería muy apreciada. Gracias !