Experiments in discrete stochastic simulation using simplified predator-prey in The Azimuth Project

# Experiments in discrete stochastic simulation using simplified predator-prey #

* toc
{: toc}

## Idea ##

Working example for author developing a C++ discrete time simulation framework.

**NOTE**: these are messing-about experiments so the graphs below don't have axis ticks labelled correctly.

## Model ##

Incredibly simplified model where both rabbits and foxes live at most 2 years, and essentially reproduce when they turn 1 year old. They hvae fertility rates specified per individual, with no tracking of males/females, and the number of babies produced by the population in total is taken as uniform rather than a more appropriate distribution. It also uses the oversimplified rule that if a fox needs to eat $c$ rabbits to survive for 1 year, then they "telepathically" arrange so that a fox eats either $c$ rabbits or no rabbits at all.

quantity | symbol
:-----------|:---------
no of foxes born this year | $f_0$
no of 1 year old foxes | $f_1$
_half_ of maximum offspring from pair of foxes | $f_f$
no of rabbits a fox needs to eat in a year to survive | $c$
no of rabbits born this year | $r_0$
no of 1 year old rabbits | $r_1$
_half_ of maximum offspring from pair of rabbits | $r_f$
maximum no rabbits vegetation can support | $K_r$

### System equations ###

The populations evolve from year to year with some stochastic equations (in discrete time):
\[\begin{aligned}
e &= clampAbove(c (f_0^t + f_1^t),r_0^t + r_1^t)\\
r_1^{t+1} &= v \quad when \quad v \geq 2 \quad where \quad v=r_0^t - clampAbove (e - r_1^t ,0)\\
f_1^{t+1} &= v \quad when \quad v \geq 2 \quad where \quad v=clampAbove (e/c,f_0^t)\\
r_0^{t+1} &= clampAbove(r_1^{t+1}  r_f  U, K_r - r_1^{t+1})\\
f_0^{t+1} &= f_1^{t+1} f_f U
\end{aligned}\]
where $U$ denotes a _fresh_ random variable uniformly distributed on $[0,1)$
for each occurrence.

## Simulation results ##

Running the system for $2^22$ different random simulation runs for various values of fox and rabbit fertility (for fixed $c$, $capR$, etc), the "fox population has not died out" probabilities after 50 years can be visualised in the below plot (where horizontal axis left to right is increasing fox fertility and vertical axis top to bottom is increasing rabbit fertility):

<div>
<img src="http://www.azimuthproject.org/azimuth/files/surivalProb.png" width="600">
</div>

The "reasonably favourable values of fertility" for foxes lie in a range of about $[2,2.6]$, and providing it's above a minimum value rabbit fertility doesn't matter.

The surival probabilities taken as a series of slices parallel to the "rabbit fertility axis":

<div>
<img src="http://www.azimuthproject.org/azimuth/files/surivalProbRabbitFertility.png" width="600">
</div>

At least part of the reason why the curves go to a horizontal line is that, with a fixed "rabbit carrying capacity", above a certain rabbit fertility level it's carrying capacity rather than rabbit fertility that determines the number of rabbits.

The surival probabilities taken as a series of slices parallel to the "fox fertility axis":

<div>
<img src="http://www.azimuthproject.org/azimuth/files/surivalProbFoxFertility.png" width="600">
</div>

This looks like it might be a gamma distribution (speculation: this might possibly be because the "surviving trajectories" are ones that don't hit any of the clamping terms, so somehow in the visualised region it's a "nice" section which is the straightforward product of uniform random variates, which may have some nice closed-form?). Running a simulation with only $2^16$ different random simulation runs has "kinks" which one couldn't tell if are genuinely significant parts of the distribution or are "under-sampling" artifacts:

<div>
<img src="http://www.azimuthproject.org/azimuth/files/surivalProbFoxFertility_2_to_16.png" width="600">
</div>

(The horizontal axis tick labels are correct here: the absolute probability of surviving 50 years is indeed below $0.006$.)

It's also interesting to see how the survival distribution evolves over time. Since there is an "absorbing barrier" at 0 (i.e., extinction) the survival probability can only decrease over time, so the distributions can be plotted on the same graph, with timesteps corresponding to consecutive curves moving downwards. Due to the drop off in scale making later curves difficult to see, these are plotted (with y-axis being probability, x-axis related to fox fertility) in blocks of 10 consecutive timesteps:

<div>
<img src="http://www.azimuthproject.org/azimuth/files/survivialFoxFert0.png" width="600">
</div>

Steps 0--9.
 
<div>
<img src="http://www.azimuthproject.org/azimuth/files/survivialFoxFert1.png" width="600">
</div>

Steps 10--19.
 
<div>
<img src="http://www.azimuthproject.org/azimuth/files/survivialFoxFert2.png" width="600">
</div>

Steps 20--29.
 
<div>
<img src="http://www.azimuthproject.org/azimuth/files/survivialFoxFert3.png" width="600">
</div>

Steps 30--39.

As can be seen, the evolution towards the skewed beta-ish distribution after 50 steps starts with some distinctly different curve shapes.

For completeness, here are the other parameters used

Fixed system parameters | Value
:-------------------------|:-------
$c$ | 200
initial $f_0$ | 100
initial $f_1$ | 0
initial $r_0$ | $2 c f_0$
initial $r_1$ | 0
$K_r$ | $5 c f_0$

## Notes ##

$2^22$ evaluations of 50 timesteps of this incredibly simple system for a $32 \times 32$ grid of system parameter values took about 3 hours on my PC.

category: software