PyMC3의 단순 동역학 모델

Question

두 개의 파라미터를 추론하기 위해 PyMC3에 동적 시스템 모델을 구성하려고합니다. 이 모델은 통상적으로 역학에 사용되는 기본 SIR이다 :

DS/DT = - R0 *의 g의 *의 S의 * I

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

어떤 도움을 주시면 감사하겠습니다. 감사 !

Answer 1

새로운 배포판을 정의하려고합니다. 다음과 같은 것. 그러나 이것은 제대로 작동하지 않으며, 내가 잘못한 것을 확실히 알지 못합니다.

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