For the sake of this question, I am attempting to port a pymc2 example from this blog post, the bit of interest being:
## SI model
from pymc import *
from numpy import *
#observed data
T = 10
susceptible_data = array([999,997,996,994,993,992,990,989,986,984])
infected_data = array([1,2,5,6,7,8,9,11,13,15])
# stochastic priors
beta = Uniform('beta', 0., 40., value=1.)
gamma = Uniform('gamma', 0., 1., value=.001)
SI_0 = Uninformative('SI_0', value=[999., 1.])
# deterministic compartmental model
@deterministic
def SI(SI_0=SI_0, beta=beta, gamma=gamma):
S = zeros(T)
I = zeros(T)
S[0] = SI_0[0]
I[0] = SI_0[1]
for i in range(1,T):
S[i] = S[i-1] - 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1])
I[i] = max(0., I[i-1] + 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1]) - gamma*I[i-1])
return S, I
S = Lambda('S', lambda SI=SI: SI[0])
I = Lambda('I', lambda SI=SI: SI[1])
# data likelihood
A = Poisson('A', mu=S, value=susceptible_data, observed=True)
B = Poisson('B', mu=I, value=infected_data, observed=True)
The part that I am stuck with is that the author has used a plain-old (backwards) difference method with a fixed step size of 0.05 to simulate S and I at each step, very sensible.
However, I cannot figure out a way to do this in pymc3 since said simulation requires running the (backwards) difference approximation over the array, which is not something that I can seem to figure out how to do in pymc3.
For example, my first approximation was to do:
import numpy as np
import pymc3 as pm
import theano
t = 10
susceptible_data = array([999,997,996,994,993,992,990,989,986,984])
infected_data = array([1,2,5,6,7,8,9,11,13,15])
def si_bdiff(si_0, beta, gamma):
s = np.zeros(t, dtype=theano.config.floatX)
f = np.zeros(t, dtype=theano.config.floatX)
s[0] = si_0[0]
f[0] = si_0[1]
for i in range(1,t):
s[i] = s[i-1] - 0.05 * beta * s[i-1] * f[i-1] / (s[i-1] + f[i-1])
f[i] = max(0., f[i-1] + 0.05 * beta * s[i-1] * f[i-1] / (s[i-1] + f[i-1]) - gamma * f[i-1])
return s, f
with pm.Model() as model:
beta = pm.Uniform('beta', lower=0, upper=40)
gamma = pm.Uniform('gamma', lower=0, upper=1)
si_0 = pm.Flat('SI_0', shape=2)
s, f = pm.Deterministic('si', si_bdiff(si_0, beta, gamma))
# data likelihood
A = pm.Poisson('A', mu=s, observed=susceptible_data)
B = pm.Poisson('B', mu=f, observed=infected_data)
trace1 = pm.sample(100000, tune=50000, cores=4, init='advi')
which, of course, failed miserably. I have been fiddling around with a number of aspects with this problem, but the above version holds my key problem, in that I don't know how to use pm.Deterministic in this manner, specifically to be able to iterate over (something?) and return what is basically a pair of arrays.
Thanks!
EDIT
This is driving me crazy, I'm so close but...
Here is the guts of my current code:
import numpy as np
import pymc3 as pm
import theano
import theano.tensor as tt
import theano.tests.unittest_tools
from scipy.optimize import fsolve
def dSdt(Scur, Icur, beta):
N = Icur + Scur
eqn = -beta * Scur * Icur / N
return eqn
def dIdt(Scur, Icur, beta, gamma):
N = Icur + Scur
eqn = beta * Scur * Icur / N - gamma * Icur
return eqn
def bkwd_euler(cur, prev, beta, gamma, h):
Scur, Icur = cur
Sprev, Iprev = prev
N = Scur + Icur
eqn0 = Scur - Sprev - h * dSdt(Scur, Icur, beta)
eqn1 = Icur - Iprev - h * dIdt(Scur, Icur, beta, gamma)
return eqn0, eqn1
T = infected_data.shape[0]
class SimulateDE(tt.Op):
itypes = [tt.dscalar, tt.dscalar, tt.dscalar, tt.dscalar]
otypes = [tt.dmatrix]
tsteps = range(T)
def perform(self, node, inputs, output):
s0, i0, beta, gamma, = inputs[0], inputs[1], inputs[2], inputs[3]
for i in self.tsteps:
if i == 0:
si = np.array((s0, i0)).reshape(1, -1)
else:
prev = si[-1]
nxt = fsolve(bkwd_euler, (100, 100), (prev, beta, gamma, 1/7))
si = np.vstack((si, nxt))
output[0][0] = si
def grad(self, inputs, g):
# this is just a placeholder, because I cannot
# figure out how grad() should work in this case
return [g[0][0], g[0][1], g[0][2], g[0][3]]
The issue I'm stuck with now is that I simply cannot figure out what grad() should be in this case, as by my thinking I should simply be able to pass g through as I have; gibberish result, but I cannot see why it shouldn't run. But alas, it does not:
theano.config.compute_test_value = 'ignore'
theano.tests.unittest_tools.verify_grad(SimulateDE(), [np.array(100.), np.array(1.), np.array(10.), np.array(1.)])
<snip>
NotImplementedError: ('input nd\nApply node that caused the error: InplaceDimShuffle{x,0}(<__main__.SimulateDE object at 0x7fdcec9594d0>.0)\nToposort index: 1\nInputs types: [TensorType(float64, vector)]\nInputs shapes: [(45, 2)]\nInputs strides: [(16, 8)]\nInputs values: [\'not shown\']\nOutputs clients: [[Elemwise{mul,no_inplace}(random_projection, InplaceDimShuffle{x,0}.0)]]\n\nBacktrace when the node is created(use Theano flag traceback.limit=N to make it longer):\n File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 2886, in _run_cell\n return runner(coro)\n File "/usr/local/lib/python3.7/dist-packages/IPython/core/async_helpers.py", line 68, in _pseudo_sync_runner\n coro.send(None)\n File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 3063, in run_cell_async\n interactivity=interactivity, compiler=compiler, result=result)\n File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 3254, in run_ast_nodes\n if (await self.run_code(code, result, async_=asy)):\n File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 3331, in run_code\n exec(code_obj, self.user_global_ns, self.user_ns)\n File "<ipython-input-108-a997c24f1e9a>", line 2, in <module>\n theano.tests.unittest_tools.verify_grad(SimulateDE(), [np.array(100.), np.array(1.), np.array(10.), np.array(1.)])\n File "/usr/local/lib/python3.7/dist-packages/theano/tests/unittest_tools.py", line 92, in verify_grad\n T.verify_grad(op, pt, n_tests, rng, *args, **kwargs)\n File "/usr/local/lib/python3.7/dist-packages/theano/gradient.py", line 1759, in verify_grad\n cost = theano.tensor.sum(t_r * o_output)\n\nHINT: Use the Theano flag \'exception_verbosity=high\' for a debugprint and storage map footprint of this apply node.', '\nThe error happened with the following inputs:', [array(100.), array(1.), array(10.), array(1.)], '\nThe value of eps is:', None, '\nThe out_type is:', None)
So that grad() seems to be getting upset about the fact that the output of bkwd_euler is a (45, 2) array, maybe? From what I can understand from the docs, the output of grad() should be a 4-dimensional vector since 4 scalar inputs, so not sure what the problem is at this point.
Any help still greatly appreciated!