Dynamical model in pymc3

Questions
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:

4

Hey if I try to go to trace = pm.sample (niter = 100, chains = 4), ti gives me error “bad intial energy”, do you know what that could have mean?

5

Bad Initial Energy is one of the worst general errors you can get, and it can mean many things. I would try running the model on a very small subset, to see if it works, and if it does, but starts failing again when the sample gets big enough, that just means Pymc3 needs those ODE solvers to chug through this model first :yum: .

1 Like
6

Now that @dpananos’s GSoC PR is almost ready: https://github.com/pymc-devs/pymc3/pull/3578, it’s time again to revive this thread :wink:

7

Remind me to post an answer once the PR is merged.