Simple Single Variable Lotka-Volterra Model in PyMC3

Questions
1

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)

2

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.