Я пытаюсь собрать модель динамической системы в PyMC3, чтобы вывести два параметра. Модель является основным SIR, обычно используемым в эпидемиологии:

дс / дт = - г0 * г * с * я

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

где r0 и g - параметры, которые будут выведены. Пока что яЯ не могу получить очень далеко вообще. Единственные примеры, которые ямы видели, как такая цепочка Маркова собирает ошибки, что приводит к тому, что рекурсия слишком глубока. Вот'мой пример кода.

# 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())

Любая помощь приветствуется. Спасибо !