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!