PyMC3中的简单动力学模型 [英] Simple Dynamical Model in PyMC3

问题描述

我正在尝试在PyMC3中建立一个动力学系统的模型，以推断出两个参数.该模型是流行病学中常用的基本SIR:

I'm trying to put together a model of a dynamical system in PyMC3, to infer two parameters. The model is the basic SIR, commonly used in epidemiology :

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

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

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

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

其中r0和g是要推断的参数.到目前为止，我根本无法做到.我所见过的将这样的马尔可夫链组合在一起的唯一示例会产生有关递归太深的错误.这是我的示例代码.

where r0 and g are parameters to be inferred. So far, I'm unable to get very far at all. The only examples I've seen of putting together a Markov chain like this yields errors about recursion being too deep. Here's my example code.

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

任何帮助将不胜感激.谢谢！

Any help would be much appreciated. Thanks !

推荐答案

我会尝试定义一个新的发行版.类似于以下内容.但是，这不是很有效，而且我不确定自己做错了什么.

I would try defining a new distribution. Something like the following. However, this is not quite working, and I'm not quite sure what I did wrong.

class SIR(Distribution): 
def __init__(self, gamma, r0,dt, std): 
    self.gamma = gamma
    self.r0 = r0
    self.std = std
    self.dt = dt

def logp(self, SI):
    r0 = self.r0 
    std = self.std 
    gamma = self.gamma 
    dt = self.dt

    S=SI[:,0]
    I=SI[:,1]

    Si = S[1:]
    Si_m1 = S[:-1]
    Ii = I[1:]
    Ii_m1 = I[:-1]

    Sdelta = (Si - Si_m1)
    Idelta = (Ii - Ii_m1)

    Sexpected_delta = dt* (-r0 * gamma * Si_m1 * Ii_m1)
    Iexpected_delta = dt * gamma * Ii_m1 *( r0 * Si_m1 - 1 )


    return (Normal.dist(Sexpected_delta, sd=std).logp(Sdelta) +
            Normal.dist(Iexpected_delta, sd=std).logp(Idelta))


with Model() as model: 
    r0 = pymc.Normal("r0", 15, sd=10)
    gamma = pymc.Normal("gamma", 0.3, 1.)
    std = .5
    dt = t[1]-t[0]


    SI = SIR('SI', gamma, r0, std,dt, observed=solution[:,:2])

    #start = pymc.find_MAP(start={'gamma' : .45, 'r0' : 17})
    trace = pymc.sample(2000, pymc.NUTS())

这篇关于PyMC3中的简单动力学模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！