Consider the following simulation problem:
Process 1: Blue-bellied mudskippers reproduce in discrete non-overlapping generations.Their population dynamics are simple; if the population size at time t is n, then the expected population size at time t+1 will be . Plot the expected population size over time.
First, we write a function corresponding that gives the next year's population size.
For example, if we were to start with 30 mudskippers, in the next generation we would have an expected population size of
We could write a loop to record the population size over a series of years:
However, notice what is actually happening mathematically. We are iterating the function p, starting from a value of 30. A preferable programming solution is to write a functional program that reflects the iterative nature of the process directly, using NestList, which works as follows.
NestList[f,a,length] creates the list
While this process is deterministic, we can also iterate stochastic functions. This is illustrated in the next example.
Process 2: Red-bellied mudskippers also reproduce in discrete non-overlapping generations,but their population dynamics depend on the weather conditions. The weather conditions can be summarized by two stochastic variables, a, uniformly distributed between 200 and 1000, and b, uniformly distributed between 0 and 10. Given the weather conditions, if the population size at time t is n, then the expected population size at time t+1 will be . Plot the expected population size over time.
There are a number of way we could code this. One would be a recursive
Do loop, as above, but that is perhaps the least elegant. Somewhat more elegant is to redefine the function p to incorporate the random variables representing the weather conditions. Below, I create a Module to do this. The
Module command can be used to define a multistep function that uses local variables.
Earlier in the tutorial, I suggested the general strategy of storing all of the simulation information. In the present example, this might be useful for a number of reasons. For example, how would the dynamics have looked - given this particular set of weather patterns - if the initial population size had been 10 instead of 30? Would it matter, in the long run, or will the population size eventually be the same regardless of where the population started?
Given the code above, there is no way to answer this question. If we try to re-run the program using a starting population size, a new series of random weather events (a and b values) will be generated. Indeed in the program above, even though we have stored all of the population size data points, we have lost the record of the random events themselves. How can we keep a record of these events, and still write elegant functional code?
We start by following the general strategy above, generating data and then applying a function. We think of the simulation itself is simply the process of generating all the necessary random events; in this case, generating a list of the random weather conditions at each time step
. We then use this list to compute the dynamics, with the
FoldList takes three arguments: a function
f, a starting value
s, and a list of parameters for each time-step,
FoldList function returns a list of successive iterations, taking the parameters as inputs:
Now, we can compare the population dynamics for different starting population sizes, given the same set of weather conditions.
It looks like the initial population size doesn't matter. Regardless of the starting size, the populations always seem to settle into the same trajectory.