<jupyter_start><jupyter_text>Lotka-Volterra toy modelThe [LotkaVolterraModel](http://pints.readthedocs.io/en/latest/toy/lotka_volterra_model.html) 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](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations).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.<jupyter_code>import pints
import pints.toy

model = pints.toy.LotkaVolterraModel()

print('Outputs: ' + str(model.n_outputs()))
print('Parameters: ' + str(model.n_parameters()))<jupyter_output>Outputs: 2
Parameters: 4<jupyter_text>The model comes pre-packaged with lynx-hare pelt count data collected by the Hudson'sBay Company in Canada in the early twentieth century, which is takenfrom [1]. The data given here corresponds to annual observationstaken from 1900-1920 (inclusive). We now plot this data.[1] <NAME>. (2009). Modeling basics. Lecture Notes for Math 442,           Texas A&M University<jupyter_code>times = model.suggested_times()
values = model.suggested_values()

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('Years since 1920')
plt.ylabel('Population size')
plt.plot(times, values)
plt.legend(['hare', 'lynx'])
plt.show()<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>plt.figure()
plt.xlim(0, 80)
plt.ylim(0, 80)
plt.xlabel('hare')
plt.ylabel('lynx')
plt.plot(values[:, 0], values[:, 1])
plt.show()<jupyter_output><empty_output><jupyter_text>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<jupyter_code>import numpy as np

# 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)

print('Running')
timer = pints.Timer()
chains = mcmc.run()
time = timer.time()
print('Done!')
import pints.plot
pints.plot.trace(chains)
plt.show()
results = pints.MCMCSummary(chains=chains, parameter_names=["a", "b", "c", "d", "sigma_1", "sigma_2"], time=time)
print(results)<jupyter_output>param    mean    std.    2.5%    25%    50%    75%    97.5%    rhat    ess    ess per sec.
-------  ------  ------  ------  -----  -----  -----  -------  ------  -----  --------------
a        0.32    0.05    0.25    0.29   0.32   0.35   0.43     1.54    23.62  0.45
b        0.13    0.02    0.10    0.12   0.13   0.14   0.18     1.56    25.11  0.48
c        1.24    0.18    0.92    1.12   1.26   1.36   1.57     1.46    23.84  0.45
d        0.38    0.05    0.28    0.34   0.38   0.41   0.48     1.46    24.49  0.46
sigma_1  0.27    0.05    0.19    0.23   0.26   0.30   0.40     1.23    52.21  0.99
sigma_2  0.41    0.08    0.26    0.37   0.41   0.46   0.57     1.05    76.69  1.45<jupyter_text>We can also compare the predictions with these values to what we found: looks like a reasonable fit.<jupyter_code># Select first chain
chain1 = chains[0]

# Remove warm up period
chain1 = chain1[500:]

# Create fine-grained list of times for plotting
n_fine = 1000
times_fine = np.linspace(min(times), max(times), n_fine)

# Plot some predictions with these samples
num_lines = 100
hare = np.zeros((n_fine, num_lines))
lynx = np.zeros((n_fine, num_lines))
for i in range(num_lines):
    temp = np.exp(model.simulate(times=times_fine, parameters=chain1[i, :4]))
    hare[:, i] = temp[:, 0]
    lynx[:, i] = temp[:, 1]
plt.figure(figsize=(15, 7))
plt.xlabel('Years since 1900')
plt.ylabel('Populations')
plt.plot(times_fine, hare, color='blue', alpha=0.01)
plt.plot(times_fine, lynx, color='orange', alpha=0.01)
plt.plot(times, values, 'o-')
plt.show()<jupyter_output><empty_output><jupyter_text>Since this is a tricky model to fit, let's use HMC to fit the same data.<jupyter_code># 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)

print('Running')
timer = pints.Timer()
chains = mcmc.run()
time = timer.time()
print('Done!')
import pints.plot
pints.plot.trace(chains)
plt.show()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>results = pints.MCMCSummary(chains=chains, parameter_names=["a", "b", "c", "d", "sigma_1", "sigma_2"], time=time)
print(results)<jupyter_output>param    mean    std.    2.5%    25%    50%    75%    97.5%    rhat    ess    ess per sec.
-------  ------  ------  ------  -----  -----  -----  -------  ------  -----  --------------
a        0.32    0.04    0.25    0.30   0.32   0.34   0.39     1.43    36.71  0.18
b        0.13    0.02    0.09    0.12   0.13   0.14   0.16     1.39    37.80  0.19
c        1.25    0.15    0.99    1.16   1.23   1.33   1.59     1.48    36.05  0.18
d        0.38    0.04    0.30    0.35   0.37   0.40   0.48     1.51    36.63  0.18
sigma_1  0.29    0.05    0.21    0.25   0.28   0.32   0.42     1.09    64.14  0.32
sigma_2  0.45    0.06    0.35    0.40   0.45   0.49   0.59     1.22    47.76  0.24
{'fraction_non_alphanumeric': 0.14222484508704633, 'fraction_numerical': 0.08070227205665388, 'mean_word_length': 3.0447494033412887, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '382050', 'license': 'RSA-MD', 'n_tokens_mistral': 2820, 'n_tokens_neox': 2499, 'n_words': 694, 'path': 'examples/toy/model-lotka-volterra.ipynb', 'repo_name': 'lisaplag/pints', 'size': '6781'}