# Simple Single Variable Lotka-Volterra Model in PyMC3

**URL:** <https://discourse.pymc.io/t/simple-single-variable-lotka-volterra-model-in-pymc3/4330>\
**Category:** Questions\
**Created:** [January 9, 2020, 3:08pm UTC](https://discourse.pymc.io/t/simple-single-variable-lotka-volterra-model-in-pymc3/4330 "2020-01-09T15:08:10Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Debbie](https://avatars.discourse-cdn.com/v4/letter/d/90ced4/32.png) [@Debbie](https://discourse.pymc.io/u/Debbie)\
**Post date:** [January 9, 2020, 3:08pm UTC](https://discourse.pymc.io/t/simple-single-variable-lotka-volterra-model-in-pymc3/4330/1 "2020-01-09T15:08:10Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [January 9, 2020, 7:38pm UTC](https://discourse.pymc.io/t/simple-single-variable-lotka-volterra-model-in-pymc3/4330/2 "2020-01-09T19:38:15Z")

</div>

Any reason not using the ODE module to model these kind of problem  
Also, you should not use find\_MAP + Metropolis: it is not a good inference approach and you likely dont get converge result.