The LotkaVolterraModel describes the relationship between two interacting species, where one preys on the other. A good description of its history and interpretation can be found on Wikipedia.
The model has 2 states $x$ and $y$, where $x$ represents a population of prey, and $y$ represents a population of predators. It is described by the ODEs:
$$ \frac{dx}{dt} = ax - bxy $$and
$$ \frac{dy}{dt} = -cy + dxy $$where $a, b, c$, and $d$ are the four model parameters.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pints
import pints.toy
import time
import pints.plot
from scipy.interpolate import interp1d
model = pints.toy.LotkaVolterraModel()
print('Outputs: ' + str(model.n_outputs()))
print('Parameters: ' + str(model.n_parameters()))
The model comes pre-packaged with lynx-hare pelt count data collected by the Hudson's Bay Company in Canada in the early twentieth century, which is taken from [1]. The data given here corresponds to annual observations taken from 1900-1920 (inclusive). We now plot this data.
[1] Howard, P. (2009). Modeling basics. Lecture Notes for Math 442, Texas A&M University
In [2]:
times = model.suggested_times() + 1900
values = model.suggested_values()
plt.figure()
plt.xlabel('Time')
plt.ylabel('Population size')
plt.plot(times, values)
plt.legend(['hare', 'lynx'])
plt.show()
In this set-up, the first state represents the prey, and the second the predators. When there is no prey, the predators begin to die out, which allows the prey population to recover.
To show the cyclical nature more clearly, these two populations are often plotted against each other:
In [3]:
plt.figure()
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.xlabel('hare')
plt.ylabel('lynx')
plt.plot(values[:, 0], values[:, 1])
plt.show()
We now use PINTS to fit the Lotka-Volterra model (with fixed initial conditions $[x,y]=[30, 4]$) to the pelts data on lynxs and hares. Previous work has showed that a multiplicative noise model is more appropriate to fit these data [2]. As such, we fit the model to the log of the series.
[2] Predator-Prey Population Dynamics: the Lotka-Volterra model in Stan. Carpenter, B. https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html
In [52]:
# Create an object with links to the model and time series
problem = pints.MultiOutputProblem(model, times, np.log(values))
# Create a log posterior
log_prior_theta = pints.UniformLogPrior(lower_or_boundaries=0, upper=2)
log_prior_sigma = pints.GaussianLogPrior(mean=0, sd=3)
log_prior = pints.ComposedLogPrior(log_prior_theta, log_prior_theta, log_prior_theta, log_prior_theta,
log_prior_sigma, log_prior_sigma)
log_likelihood = pints.GaussianLogLikelihood(problem)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Run MCMC on the noisy data
x0 = [[0.43, 0.16, 0.9, 0.27, 0.28, 0.26]] * 4
mcmc = pints.MCMCController(log_posterior, 4, x0)
mcmc.set_max_iterations(4000)
#mcmc.set_log_to_screen(False)
start = time.time()
print('Running')
chains = mcmc.run()
print('Done!')
end = time.time()
diff=end-start
In [53]:
import pints.plot
pints.plot.trace(chains)
plt.show()
In [54]:
results = pints.MCMCSummary(chains=chains, parameter_names=["a", "b", "c", "d", "sigma_1", "sigma_2"], time=diff)
print(results)
We can also compare the predictions with these values to what we found: looks like a reasonable fit.
In [45]:
# Select first chain
chain1 = chains[0]
# Remove burn-in
chain1 = chain1[500:]
# Plot some predictions with these samples
num_lines = 100
hare = np.zeros((len(times), num_lines))
lynx = np.zeros((len(times), num_lines))
for i in range(num_lines):
temp = np.exp(model.simulate(times=times, parameters=chain1[i, :4]))
hare[:, i] = temp[:, 0]
lynx[:, i] = temp[:, 1]
plt.plot(times, hare, color="b", alpha=0.1)
plt.plot(times, lynx, color="orange", alpha=0.1)
plt.plot(times, values, 'o-')
plt.show()
Since this is a tricky model to fit, let's use HMC to fit the same data.
In [46]:
# Run MCMC on the noisy data
x0 = [[0.43, 0.16, 0.9, 0.27, 0.28, 0.26]] * 4
mcmc = pints.MCMCController(log_posterior, 4, x0, method=pints.HamiltonianMCMC)
mcmc.set_max_iterations(200)
mcmc.set_log_interval(1)
for sampler in mcmc.samplers():
sampler.set_leapfrog_step_size([0.1, 0.01, 0.1, 0.03, 0.05, 0.05])
sampler.set_leapfrog_steps(10)
start = time.time()
print('Running')
chains = mcmc.run()
print('Done!')
end = time.time()
diff=end-start
In [47]:
import pints.plot
pints.plot.trace(chains)
plt.show()
We get similar results as with adaptive covariance; except, the efficiency suffers due to having to calculate the sensitivities. Overall, for this problem, adaptive covariance performs favourably.
In [51]:
results = pints.MCMCSummary(chains=chains, parameter_names=["a", "b", "c", "d", "sigma_1", "sigma_2"], time=diff)
print(results)