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.
import numpy as np
import matplotlib.pyplot as plt
import pints
import pints.toy
model = pints.toy.LotkaVolterraModel()
print('Outputs: ' + str(model.n_outputs()))
print('Parameters: ' + str(model.n_parameters()))
Outputs: 2 Parameters: 4
We can get some suggested parameters and plot some data!
times = model.suggested_times()
parameters = model.suggested_parameters()
print(parameters)
plt.figure()
plt.xlabel('Time')
plt.ylabel('Population size')
plt.plot(times, model.simulate(parameters, times))
plt.legend(['prey', 'predators'])
plt.show()
[3 2 3 2]
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:
values = model.simulate(parameters, times)
plt.figure()
plt.xlim(0, 10)
plt.ylim(0, 10)
plt.xlabel('Prey')
plt.ylabel('Predators')
plt.plot(values[:, 0], values[:, 1])
plt.show()
We can explore the model further by varying the initial conditions:
plt.figure()
plt.xlim(0, 8)
plt.ylim(0, 8)
plt.xlabel('Prey')
plt.ylabel('Predators')
for x in [1.5, 1.6, 2, 2.5, 3, 4, 5]:
model.set_initial_conditions([x, x])
values = model.simulate(parameters, times)
plt.plot(values[:, 0], values[:, 1])
plt.show()
To see how sensitive the system is to different parameters, we can set up and run an MCMC routine:
model.set_initial_conditions([2, 2])
values = model.simulate(parameters, times)
sigma = 0.05
noisy_values = values + np.random.normal(0, sigma, values.shape)
plt.figure()
plt.xlabel('Time')
plt.ylabel('Population size')
plt.plot(times, noisy_values)
plt.plot(times, values)
plt.show()
# Create an object with links to the model and time series
problem = pints.MultiOutputProblem(model, times, noisy_values)
# Create a log posterior
log_prior = pints.UniformLogPrior([1, 1, 1, 1], [6, 6, 6, 6])
log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Run MCMC on the noisy data
x0 = [[4, 1, 2, 3]]*3
mcmc = pints.MCMCController(log_posterior, 3, x0)
mcmc.set_max_iterations(4000)
#mcmc.set_log_to_screen(False)
print('Running')
chains = mcmc.run()
print('Done!')
Running Using Adaptive covariance MCMC Generating 3 chains. Running in sequential mode. Iter. Eval. Accept. Accept. Accept. Time m:s 0 3 0 0 0 0:00.0 1 6 0 0.5 0 0:00.0 2 9 0 0.667 0 0:00.0 3 12 0.25 0.5 0 0:00.0 20 63 0.381 0.524 0.286 0:00.1 40 123 0.439 0.415 0.415 0:00.2 60 183 0.377 0.344 0.377 0:00.2 80 243 0.321 0.333 0.296 0:00.3 100 303 0.257 0.297 0.238 0:00.3 120 363 0.214876 0.273 0.198 0:00.4 140 423 0.184 0.234 0.17 0:00.5 160 483 0.161 0.205 0.149 0:00.5 180 543 0.144 0.188 0.133 0:00.6 Initial phase completed. 200 603 0.129 0.174 0.119403 0:00.7 220 663 0.127 0.172 0.122 0:00.7 240 723 0.133 0.191 0.129 0:00.8 260 783 0.176 0.192 0.161 0:00.9 280 843 0.196 0.199 0.185 0:00.9 300 903 0.216 0.229 0.203 0:01.0 320 963 0.227 0.252 0.218 0:01.1 340 1023 0.246 0.264 0.228739 0:01.1 360 1083 0.252 0.258 0.241 0:01.2 380 1143 0.255 0.268 0.244 0:01.3 400 1203 0.257 0.267 0.254 0:01.3 420 1263 0.261 0.263658 0.259 0:01.4 440 1323 0.263 0.260771 0.254 0:01.4 460 1383 0.271 0.262 0.249 0:01.5 480 1443 0.262 0.258 0.247 0:01.6 500 1503 0.257485 0.265 0.243513 0:01.6 520 1563 0.255 0.263 0.251 0:01.7 540 1623 0.253 0.266 0.253 0:01.7 560 1683 0.254902 0.267 0.253 0:01.8 580 1743 0.258 0.272 0.251 0:01.9 600 1803 0.256 0.263 0.245 0:01.9 620 1863 0.252818 0.259 0.245 0:02.0 640 1923 0.253 0.254 0.245 0:02.1 660 1983 0.251 0.248 0.245 0:02.1 680 2043 0.253 0.248 0.244 0:02.2 700 2103 0.251 0.243 0.243 0:02.2 720 2163 0.247 0.245 0.239 0:02.3 740 2223 0.245614 0.244 0.235 0:02.4 760 2283 0.243 0.244 0.23 0:02.4 780 2343 0.242 0.239 0.224 0:02.5 800 2403 0.241 0.24 0.22 0:02.6 820 2463 0.242 0.24 0.218 0:02.6 840 2523 0.237 0.235434 0.21522 0:02.7 860 2583 0.233 0.233 0.213 0:02.8 880 2643 0.238 0.238 0.211 0:02.9 900 2703 0.235 0.238 0.211 0:02.9 920 2763 0.236 0.24 0.212 0:03.0 940 2823 0.238 0.24017 0.214 0:03.1 960 2883 0.234 0.24 0.219 0:03.1 980 2943 0.234 0.243629 0.218 0:03.2 1000 3003 0.231 0.243 0.219 0:03.3 1020 3063 0.233 0.24 0.224 0:03.4 1040 3123 0.233 0.24 0.221902 0:03.4 1060 3183 0.236 0.24 0.221 0:03.5 1080 3243 0.237 0.239593 0.221 0:03.6 1100 3303 0.238 0.242 0.222 0:03.6 1120 3363 0.238 0.243 0.224 0:03.7 1140 3423 0.237511 0.245 0.223 0:03.7 1160 3483 0.235 0.246 0.221 0:03.8 1180 3543 0.233 0.249 0.223 0:03.9 1200 3603 0.234 0.246 0.22398 0:03.9 1220 3663 0.237 0.244 0.221 0:04.1 1240 3723 0.233 0.242 0.218 0:04.1 1260 3783 0.234 0.239 0.219 0:04.2 1280 3843 0.233 0.24 0.219 0:04.3 1300 3903 0.23 0.239 0.223 0:04.3 1320 3963 0.23 0.238 0.221 0:04.4 1340 4023 0.233 0.239 0.223 0:04.5 1360 4083 0.231 0.238795 0.226 0:04.5 1380 4143 0.232 0.237 0.226 0:04.6 1400 4203 0.23 0.238 0.226 0:04.7 1420 4263 0.229 0.24 0.228 0:04.7 1440 4323 0.232 0.243 0.229 0:04.8 1460 4383 0.231 0.247091 0.229 0:04.9 1480 4443 0.231 0.246 0.229 0:05.0 1500 4503 0.232 0.245 0.229 0:05.1 1520 4563 0.229 0.246 0.227 0:05.2 1540 4623 0.229721 0.245 0.227 0:05.2 1560 4683 0.233 0.243 0.228 0:05.3 1580 4743 0.231 0.242 0.227704 0:05.4 1600 4803 0.232 0.242 0.229 0:05.6 1620 4863 0.231 0.244 0.23 0:05.6 1640 4923 0.233 0.245582 0.23 0:05.7 1660 4983 0.232 0.245 0.230584 0:05.8 1680 5043 0.233 0.244 0.233 0:05.9 1700 5103 0.236 0.242 0.233 0:06.0 1720 5163 0.236 0.242882 0.233 0:06.1 1740 5223 0.236 0.243 0.231 0:06.2 1760 5283 0.236 0.241 0.229983 0:06.3 1780 5343 0.236 0.242 0.232 0:06.4 1800 5403 0.234 0.242643 0.232 0:06.4 1820 5463 0.234 0.243 0.231 0:06.6 1840 5523 0.234 0.242 0.232 0:06.7 1860 5583 0.232 0.241 0.232 0:06.8 1880 5643 0.23 0.242 0.232 0:06.8 1900 5703 0.231 0.243 0.232 0:06.9 1920 5763 0.231 0.242582 0.233 0:07.0 1940 5823 0.229 0.244204 0.235 0:07.1 1960 5883 0.229 0.243 0.234 0:07.1 1980 5943 0.229 0.242 0.232 0:07.2 2000 6003 0.228 0.242 0.231 0:07.3 2020 6063 0.228 0.241 0.232 0:07.4 2040 6123 0.227 0.243 0.23 0:07.4 2060 6183 0.225 0.242 0.23 0:07.5 2080 6243 0.226 0.241 0.231 0:07.6 2100 6303 0.227 0.24 0.23 0:07.7 2120 6363 0.227 0.239 0.229 0:07.8 2140 6423 0.227 0.239 0.23 0:07.9 2160 6483 0.228 0.24 0.23 0:07.9 2180 6543 0.228 0.241 0.229 0:08.0 2200 6603 0.228 0.242 0.229 0:08.0 2220 6663 0.228 0.241 0.229 0:08.1 2240 6723 0.227 0.241 0.227577 0:08.2 2260 6783 0.227 0.24 0.227 0:08.3 2280 6843 0.228 0.24 0.227 0:08.4 2300 6903 0.227 0.239 0.228 0:08.5 2320 6963 0.227 0.239 0.227 0:08.5 2340 7023 0.228 0.239214 0.226399 0:08.6 2360 7083 0.226 0.239 0.226 0:08.7 2380 7143 0.226 0.24 0.227 0:08.8 2400 7203 0.226 0.2399 0.226 0:08.9 2420 7263 0.227 0.239 0.228005 0:09.0 2440 7323 0.227 0.238 0.229 0:09.1 2460 7383 0.227 0.238 0.229 0:09.2 2480 7443 0.227 0.237 0.23 0:09.2 2500 7503 0.227509 0.237 0.229 0:09.3 2520 7563 0.226 0.236 0.227 0:09.4 2540 7623 0.226 0.235734 0.227 0:09.4 2560 7683 0.225 0.235 0.227 0:09.5 2580 7743 0.224 0.235 0.227 0:09.6 2600 7803 0.225298 0.235 0.229 0:09.6 2620 7863 0.224 0.237 0.228 0:09.7 2640 7923 0.225 0.236 0.228 0:09.8 2660 7983 0.224 0.236 0.228 0:09.8 2680 8043 0.224 0.234614 0.227527 0:09.9 2700 8103 0.225472 0.234 0.226953 0:10.0 2720 8163 0.225 0.233 0.227 0:10.1 2740 8223 0.22583 0.234 0.228 0:10.2 2760 8283 0.226 0.235 0.229 0:10.3 2780 8343 0.227 0.234 0.228 0:10.3 2800 8403 0.227 0.232774 0.228 0:10.4 2820 8463 0.227 0.233 0.227 0:10.4 2840 8523 0.228 0.232 0.228 0:10.5 2860 8583 0.227 0.233 0.22964 0:10.6 2880 8643 0.228 0.233 0.22874 0:10.6 2900 8703 0.229 0.232 0.228 0:10.7 2920 8763 0.229 0.233 0.228 0:10.7 2940 8823 0.228 0.232 0.227 0:10.8 2960 8883 0.227 0.231003 0.227 0:10.9 2980 8943 0.226 0.23 0.228 0:10.9 3000 9003 0.226 0.231 0.227 0:11.0 3020 9063 0.225422 0.23 0.228 0:11.1 3040 9123 0.225 0.231 0.228 0:11.1 3060 9183 0.225 0.231 0.228 0:11.2 3080 9243 0.224927 0.231 0.228 0:11.3 3100 9303 0.224 0.231 0.227 0:11.3 3120 9363 0.224 0.232 0.227 0:11.4 3140 9423 0.224 0.231455 0.227 0:11.5 3160 9483 0.224 0.232 0.227 0:11.5 3180 9543 0.224 0.231 0.226 0:11.6 3200 9603 0.225 0.231 0.226 0:11.6 3220 9663 0.226 0.231 0.226 0:11.7 3240 9723 0.226 0.232 0.226 0:11.8 3260 9783 0.225391 0.232 0.227 0:11.9 3280 9843 0.226 0.232 0.227 0:12.0 3300 9903 0.226 0.231445 0.227 0:12.1 3320 9963 0.225 0.231 0.228 0:12.1 3340 10023 0.224 0.232 0.228 0:12.2 3360 10083 0.224 0.232 0.228801 0:12.3 3380 10143 0.224 0.232 0.229 0:12.3 3400 10203 0.224 0.233 0.229 0:12.4 3420 10263 0.224 0.234 0.228 0:12.5 3440 10323 0.224 0.233653 0.227 0:12.5 3460 10383 0.223 0.235 0.227102 0:12.6 3480 10443 0.223 0.235 0.227 0:12.6 3500 10503 0.222 0.235647 0.228 0:12.7 3520 10563 0.222096 0.235 0.228 0:12.8 3540 10623 0.223 0.235 0.227 0:12.8 3560 10683 0.223 0.233 0.227 0:12.9 3580 10743 0.223 0.233 0.227 0:12.9 3600 10803 0.224 0.234 0.227 0:13.0 3620 10863 0.223 0.235 0.228 0:13.1 3640 10923 0.223 0.235 0.229 0:13.1 3660 10983 0.224 0.235 0.23 0:13.2 3680 11043 0.223 0.236 0.23 0:13.3 3700 11103 0.222 0.235612 0.23 0:13.4 3720 11163 0.221 0.235 0.23 0:13.5 3740 11223 0.221 0.235 0.23 0:13.6 3760 11283 0.22 0.234778 0.23 0:13.6 3780 11343 0.22 0.235123 0.23 0:13.7 3800 11403 0.219 0.235 0.229 0:13.8 3820 11463 0.219 0.234 0.23 0:13.9 3840 11523 0.219 0.234 0.23 0:14.0 3860 11583 0.219 0.234 0.23 0:14.0 3880 11643 0.218 0.234 0.23 0:14.1 3900 11703 0.218 0.233 0.23 0:14.2 3920 11763 0.218 0.233 0.23 0:14.2 3940 11823 0.218 0.232 0.23 0:14.3 3960 11883 0.218 0.232517 0.229235 0:14.3 3980 11943 0.218 0.232 0.229 0:14.4 4000 12000 0.21875 0.232 0.2295 0:14.5 Halting: Maximum number of iterations (4000) reached. Done!
# Check convergence using rhat criterion
print('R-hat:')
print(pints.rhat_all_params(chains))
import pints.plot
pints.plot.trace(chains, ref_parameters=parameters)
plt.show()
R-hat: [1.078403673118607, 1.0477124994282836, 1.0651132349574168, 1.0629509070339487]
We can also compare the predictions with these values to what we found:
# Select first chain
chain1 = chains[0]
# Remove burn-in
chain1 = chain1[2500:]
# Plot some predictions with these samples
pints.plot.series(chain1, problem)
plt.show()