Dynamical model in pymc3


#1

Hi everyone,

I am trying to infer the parameters of a Lotka-Volterra equations using Pymc3. I was following this Stan example. Is there a working example of how to do similar thing in Pymc3? I can’t seem to find any. I did the following the code but it is not working

def rhs_ode(y, a, b, c, d, dt):
    yprime = T.zeros_like(y)
    yprime = T.set_subtensor(yprime[0], y[0] + dt * (a*y[0] - b*y[0]*y[1]))
    yprime = T.set_subtensor(yprime[1], y[1] + dt * (-c*y[1] + d*y[0]*y[1]))
    return yprime

#Model
LV_Model = pm.Model()
dt = t[1] - t[0]
with LV_Model:
    #priors
    alpha = pm.Normal('alpha', mu=1, sd=0.5)
    gamma = pm.Normal('gamma', mu=1, sd=0.5)
    beta = pm.Normal('beta', mu=0.05, sd=0.05)
    delta = pm.Normal('delta', mu=0.05, sd=0.05)
    sigma = pm.Lognormal('sigma', mu=-1, tau=1, shape=2)
    dt = T.fscalar('dt')
    steps = T.iscalar('steps')

    #Initial Conditions
    Z0 = pm.Lognormal('Z0', mu=np.log(10), tau=1, shape=2)

   
    # Symbolic loop through Euler updates
    xout, updates = theano.scan(fn= rhs_ode,
                            outputs_info=Z0,
                            non_sequences=[alpha, beta, gamma, delta, dt],
                            n_steps=steps)

    simulation = theano.function(inputs= [Z0, alpha, beta, gamma, delta, dt, steps],
                           outputs=xout,
                           updates=updates,
                           allow_input_downcast=True)
   
    Y_obs = pm.Lognormal('Y_obs', mu=simulation(Z0, alpha, beta, gamma, delta, dt, steps), tau=sigma,observed=data)

I get the following error

TypeError: ('Bad input argument to theano function with name "<ipython-input-30-cc4f13908e04>:48"  at index 0(0-based)', 'Expected an array-like object, but found a Variable: maybe you are trying to call a function on a (possibly shared) variable instead of a numeric array?')

The problem is, I think, with how I use theano scan and its function, and variable types. How could I make this work?


#2

This would be something difficult to do, as we dont have a ODE solver implemented yet in PyMC3. cc @michaelosthege who also works on ODE.

To write down the model in PyMC3 is possible, and you are actually quite close already:

First check the ODE function is correctly implemented in theano

theano.config.test_values = 'raise'

alpha = tt.dscalar("alpha")
beta = tt.dscalar("beta")
gamma = tt.dscalar("gamma")
delta = tt.dscalar("delta")
Z0 = tt.vector("Z0")
dt = tt.dscalar("dt")
steps = tt.iscalar("steps")

def rhs_ode(y, a, b, c, d, dt):
    yprime = tt.zeros_like(y)
    yprime = tt.set_subtensor(
        yprime[0], y[0] + dt * (a * y[0] - b * y[0] * y[1]))
    yprime = tt.set_subtensor(
        yprime[1], y[1] + dt * (-c * y[1] + d * y[0] * y[1]))
    return yprime


# Symbolic loop through Euler updates
xout, updates = theano.scan(fn=rhs_ode,
                            outputs_info=Z0,
                            non_sequences=[alpha, beta, gamma, delta, dt],
                            n_steps=steps)

simulation = theano.function(inputs=[Z0, alpha, beta, gamma, delta, dt, steps],
                                 outputs=xout,
                                 updates=updates,
                                 allow_input_downcast=True)

data = simulation(np.random.randn(2), 1, 1, .05, .05, .05, 5000)
plt.plot(data[:, 0], data[:, 1]);

Using the simulated data as above, now write down the model in PyMC3:

with pm.Model() as LV_Model:
    # priors
    alpha = pm.Normal('alpha', mu=1, sd=0.5)
    gamma = pm.Normal('gamma', mu=1, sd=0.5)
    beta = pm.Normal('beta', mu=0.05, sd=0.05)
    delta = pm.Normal('delta', mu=0.05, sd=0.05)
    sigma = pm.HalfNormal('sigma', sd=1, shape=2)

    # Initial Conditions
    Z0 = pm.Normal('Z0', mu=0, sd=1, shape=2)

    # Symbolic loop through Euler updates
    xout, updates = theano.scan(fn=rhs_ode,
                                outputs_info=Z0,
                                non_sequences=[alpha, beta, gamma, delta, .05],
                                n_steps=5000)

    Y_obs = pm.Normal('Y_obs', mu=xout, sd=sigma, observed=data)

You can now sample (slowly) from this model.
Notice that I modify the dt and steps, previously you are setting them as a theano tensor, but I doubted you are going to infer this two from the data. So I am setting it to a numpy input here instead.

Hope this helps! If you manage to port it into PyMC3 it would be a welcome contribution to the doc :wink:


AttributeError: module 'pymc3' has no attribute 'Simulator'