2
我試圖在PyMC3中組合一個動態系統模型來推斷兩個參數。該模型是基本SIR,在流行病學常用:PyMC3中的簡單動態模型
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())
任何幫助將不勝感激。謝謝 !
您是否能夠取得進展? –