python动力学建模与仿真_PyMC3中的简单动力学模型

我试图在PyMC3中建立一个动力学系统的模型，来推断两个参数。该模型是流行病学中常用的基本SIR：

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

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

任何帮助都将不胜感激。谢谢！在