# Modeling and Simulation in Python

Rabbit example

``````

In [1]:

%matplotlib inline

from modsim import *

``````

## Rabbit Redux

This notebook starts with a version of the rabbit population growth model and walks through some steps for extending it.

In the original model, we treat all rabbits as adults; that is, we assume that a rabbit is able to breed in the season after it is born. In this notebook, we extend the model to include both juvenile and adult rabbits.

As an example, let's assume that rabbits take 3 seasons to mature. We could model that process explicitly by counting the number of rabbits that are 1, 2, or 3 seasons old. As an alternative, we can model just two stages, juvenile and adult. In the simpler model, the maturation rate is 1/3 of the juveniles per season.

To implement this model, make these changes in the System object:

1. Replace `p0` with two initial populations: `juvenile_pop0` and `adult_pop0`, with values `0` and `10`.

2. Add an additional variable, `mature_rate`, with the value `0.33`.

``````

In [2]:

system = System(t0 = 0,
t_end = 30,
juvenile_pop0 = 0,
birth_rate = 0.9,
mature_rate = 0.33,
death_rate = 0.5)

system

``````
``````

Out[2]:

value

t0
0.00

t_end
30.00

juvenile_pop0
0.00

10.00

birth_rate
0.90

mature_rate
0.33

death_rate
0.50

``````

Now update `run_simulation` with the following changes:

1. Add a second TimeSeries, named `juveniles`, to keep track of the juvenile population, and initialize it with `juvenile_pop0`.

2. Inside the for loop, compute the number of juveniles that mature during each time step.

3. Also inside the for loop, add a line that stores the number of juveniles in the new `TimeSeries`. For simplicity, let's assume that only adult rabbits die.

4. During each time step, subtract the number of maturations from the juvenile population and add it to the adult population.

``````

In [3]:

def run_simulation(system):
"""Runs a proportional growth model.

Adds TimeSeries to `system` as `results`.

system: System object with t0, t_end, p0,
birth_rate and death_rate
"""
juveniles = TimeSeries()
juveniles[system.t0] = system.juvenile_pop0

for t in linrange(system.t0, system.t_end):
maturations = system.mature_rate * juveniles[t]

else:
market = 0

juveniles[t+1] = juveniles[t] + births - maturations

system.juveniles = juveniles

``````

Test your changes in `run_simulation`:

``````

In [4]:

run_simulation(system)

``````
``````

Out[4]:

value

0
10.000000

1
5.000000

2
5.470000

3
6.209900

4
7.057723

5
8.021560

6
9.117031

7
10.362107

8
11.777219

9
13.385586

10
15.213601

11
17.291261

12
19.652658

13
22.336542

14
25.386953

15
28.853947

16
32.794414

17
34.478600

18
36.487431

19
37.893339

20
39.401924

21
40.546917

22
41.694992

23
42.613800

24
43.495581

25
44.226171

26
44.907656

27
45.485241

28
46.014130

29
46.469075

30
46.880673

31
47.238170

``````

Next, update `plot_results` to plot both the adult and juvenile `TimeSeries`.

``````

In [5]:

def plot_results(system, title=None):
"""Plot the estimates and the model.

system: System object with `results`
"""
newfig()
plot(system.juveniles, 'gs-', label='juveniles')
decorate(xlabel='Season',
ylabel='Rabbit population',
title=title)

``````

And test your updated version of `plot_results`.

``````

In [6]:

plot_results(system, title='Proportional growth model')

``````
``````

``````

This notebook demonstrates the steps we recommend for starting your project:

1. Start with one of the examples from the book, either by copying a notebook or pasting code into a new notebook. Get the code working before you make any changes.

2. Make one small change, and run the code again.

3. Repeat step 2 until you have a basic implementation of your model.

If you start with working code that you understand and make small changes, you can avoid spending a lot of time debugging.

One you have a basic model working, you can think about what metrics to measure, what parameters to sweep, and how to use the model to predict, explain, or design.

### Bonus question

Suppose you only have room for 30 adult rabbits. Whenever the adult population exceeds 30, you take any excess rabbits to market (as pets for kind children, of course). Modify `run_simulation` to model this strategy. What effect does it have on the behavior of the system?

What happens if rabbits only take one season to mature? Change `mature_rate` to 1.0 and run the simulation again.

``````

In [ ]:

``````