Hi all,

I trying to create a dynamical model for my PhD. I am new to PyMC3 (and in fact, Bayesian modelling in general) so wanted to start with something simple: a Lotka-Volterra (predator-prey) model.

$$ dx/dt = \alpha x = \beta xy$$ $$dy/dt = \delta xy = \gamma y$$

I have managed to make the model run, however the predicted value fits very badly with the data provided, even with 10000 samples. I assume there is either an error in my code or in my Bayesian logic, but I can’t work out what! Here is my code:

#test value for alpha

alpha = 0.55

#Define time domain

t = np.linspace(1900,1920,21)

#Define Lotka-Volterra Equations as a function of alpha

def dy_dt(y, t, alpha):

“”“Right hand side of Lotka-Volterra equation; y = (hares, lynx)”""

beta = 0.028

delta = 0.026

gamma = 0.84

# Unpack y

h, l = y

```
# Compute derivatives
dh_dt = alpha * h - beta * l * h
dl_dt = delta * h * l - gamma * l
return np.array([dh_dt, dl_dt])
```

#Initial Condition

y0 = np.array([30, 4])

#Simulated value for Lynx and Hare populations

y_sim = odeint(dy_dt, y0, t, args=(alpha,))

# Take out Hares and Lynxes

h = y_sim[:,0]

l = y_sim[:,1]

#Data

Year = np.arange(1900,1921,1)

Lynx = np.array([4.0, 6.1, 9.8, 35.2, 59.4, 41.7, 19.0, 13.0, 8.3, 9.1, 7.4,

8.0, 12.3, 19.5, 45.7, 51.1, 29.7, 15.8, 9.7, 10.1, 8.6])

Hare = np.array([30.0, 47.2, 70.2, 77.4, 36.3, 20.6, 18.1, 21.4, 22.0, 25.4,

27.1, 40.3, 57.0, 76.6, 52.3, 19.5, 11.2, 7.6, 14.6, 16.2, 24.7])

#Define the Data Matrix

Y = np.vstack((Hare, Lynx)).T

#Create MCMC model

mcmc_model = pm.Model()

with mcmc_model:

#Priors

alpha = pm.Normal(‘alpha’, mu = 0, sigma = 100)

```
#Expected Value of Outcome
mu = y_sim
#Define the start point
start = map_estimate = pm.find_MAP()
#Define sampling method (otherwise it will just default to NUTS)
step = pm.Metropolis()
#Likelihood (sampling distribution) of observations
Y_obs = pm.Lognormal('Y_obs', mu=mu, sigma = 2, observed = Y)
#MCMC Magic!
trace = pm.sample(10000, step=step, start=start, chains = 5)
```

burnin = 100

trace = trace[burnin:]

pm.traceplot(trace)

my_summary = pm.summary(trace).round(2)

pd.DataFrame(my_summary)