In [1]:
!date
Mon Jan 27 10:19:59 PST 2014
In [2]:
%matplotlib inline
import pymc, pymc as pm, numpy as np
From: [email protected] [mailto:[email protected]] On Behalf Of Quentin CAUDRON
Sent: Monday, January 20, 2014 4:16 PM
To: [email protected]
Subject: [pymc] SIR model with PyMC3 - dynamical / compartmental

Hi,

I'm trying to put together some parameter inference for a compartmental epidemiological model. The simplest form looks like this :

\begin{align} dS/dt &= - r * g * S * I\\ dI/dt &= g * I ( r * S - 1 ) \end{align}

I'd like to use PyMC3 to infer the parameters r and g.

odeint won't accept pymc variables, so that won't work. I tried using manual integration ( just a basic forward Euler works here ), but trying to do arithmetic on pymc variables in a loop gives me ValueError: setting an array element with a sequence. :

In [3]:
t = np.linspace(0,1,20)
solution = np.zeros_like(t)
In [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]
    S = np.zeros_like(solution)
    I = np.zeros_like(solution)
    S[0] = 0.99
    I[0] = 0.01

    for i in range(len(t[1:])) :
        S[i] = S[i-1] + dt * (-r0 * gamma * S[i-1] * I[i-1])
        I[i] = I[i-1] + dt * ( r0 * gamma * S[i-1] * I[i-1] - gamma * I[i-1])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-5-c6feb0e9eed4> in <module>()
     11 
     12     for i in range(len(t[1:])) :
---> 13         S[i] = S[i-1] + dt * (-r0 * gamma * S[i-1] * I[i-1])
     14         I[i] = I[i-1] + dt * ( r0 * gamma * S[i-1] * I[i-1] - gamma * I[i-1])

ValueError: setting an array element with a sequence.

The problem is S is zeros too much like t, specifically an array of floats.

In [6]:
S.dtype
Out[6]:
dtype('float64')

Forcing S.dtype to be object sorts that out.

In [7]:
with pm.Model() as model :
    r0 = pm.Normal("r0", 15, sd=10)
    gamma = pm.Uniform("gamma", 0.3, 1.)
    
    # Use forward Euler to solve
    dt = t[1] - t[0]
    S = np.zeros_like(solution, dtype=object)  # FIX #1
    I = np.zeros_like(solution, dtype=object)  # FIX #2
    S[0] = 0.99
    I[0] = 0.01

    for i in range(len(t[1:])) :
        S[i] = S[i-1] + dt * (-r0 * gamma * S[i-1] * I[i-1])
        I[i] = I[i-1] + dt * ( r0 * gamma * S[i-1] * I[i-1] - gamma * I[i-1])

Success!