Dynamical model in pymc3

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