# monte carlo simulation using exponential distributions

I'm trying to simulate a stochastic model of deterministic exponential population growth, where $dN/dt = rN$ where $N$ is population size and $r$ is rate ($t$ time). I'm assuming there's no carrying capacity. This page (http://cnr.lwlss.net/DiscreteStochasticLogistic/) suggests this algorithm for simulating growth on an interval $[0, t_{end}]$:

1. start at $t = 0$ with initial population size
2. draw next time for birth event, $\delta t \sim Exponential(rN(1 - N/K))$ ($K$ is carrying capacity)
3. increment population size, $N = N + 1$
4. set $t = t + \delta t$
5. if $t > t_{end}$ then quit, otherwise go to step 2.

since I don't have a carrying capacity $K$, I assume it's infinite, so the next birth time is $\delta t \sim Exponential(rN)$. Is that correct?

when I run the simulation this way it doesn't at all give similar results to $N(t) = P_0 e^{rt}$ (where $P_0$ is initial population size). even averaging over many iterations it seems to give different growth curves.

below is my code and the result of the simulation. red curve is deterministic exponential growth, black curves are simulation using exponential distribution. they clearly don't match.

import numpy as np
import matplotlib.pylab as plt
def sim(rate, start, end, init):
N = 200
finalsizes = []
results = []
for n in range(N):
size = init
curr_t = 0
times = [curr_t]
sizes = [init]
new_rate = rate
while curr_t <= end:
# simulate next birth time. the scale
# parameter is inversely proportional to population
# size
new_rate = 1/float(new_rate * size)
div_time = np.random.exponential(scale=new_rate)
curr_t += div_time
if curr_t > end:
# if we exceed time interval, quit
break
times.append(curr_t)
# increase population size
size += 1
sizes.append(size)
finalsizes.append([times, sizes])
return finalsizes

# run simulation and plot results
init = 20
start = 0
end = 20
rate = 1
finalsizes = sim(rate, start, end, init)
plt.figure()
allsizes = []
for f in finalsizes:
allsizes.append(f[1][-1])
plt.plot(f[0], f[1], color="k", alpha=0.5)
times = np.arange(0, end + 1)
plt.plot(times, init*np.power(2, rate * times), color="r")
plt.xlabel("time")
plt.ylabel("size of population")
print "mean final size: ", np.mean(allsizes)
plt.show()


response to whuber's excellent answer: i don't understand why I have to specify the population size in advance. my simulation is meant to ask: in a given amount of time, what will be the variability in population size assuming exponential growth? (and not, how long would it take on average to make a population of size $N$, which is what whuber's simulation seems to be doing).

also, i don't think it's correct that i "haven't yet done enough simulations to appreciate what they are telling you". i updated my simulation to plot 1000 runs. as you can see, with my time-based stop condition, the results consistently underestimate the deterministic exponential growth population size.

my reasoning was that if i simulate $N$ runs for a duration $t$, then as $N \rightarrow \infty$, the average population size i get from my simulation should be an unbiased estimate of the population size based on the deterministic exponential growth model for $t$, i.e. $2^{rt}$. for example, if i simulate 3 time steps (starting with a single individual), i would expect the population size to sometimes be greater than 8 in the simulation and sometimes less than 8, and the average to converge to 8. it seems to be like this should be true even if i start with a very small population, as long as i simulate enough runs. is this incorrect? what is wrong with this reasoning? the simulations don't support this although i expected it to be true. it seems like there has to be a flaw in my reasoning and/or simulation.

update 2: fixed simulation where the scale parameter to exponential distribution decreases with population size (inversely proportional to it) and initial population size is 10. it still grossly underestimates exponential growth.

• Could you post some pictures of your result, alongside the proposed growth curve? – Alex R. Apr 7 '16 at 22:20
• @AlexR. posted graph and code for simulation – mvd Apr 7 '16 at 23:20
• If I understand the code correctly, it has at least two big problems. The first is the rate calculation: you seem to multiply the original rate by the initial population, then by the initial population+1, etc. With an initial population of $1$ and and initial value of $r$, this will give you a rate of $n!r$ at step $n$ rather than $nr$. The only reason you don't get a ridiculous explosion in the populations is that you appear to use the rate as a scale, which is the reciprocal of the rate. In effect, your rate is really going rapidly to zero, so the population levels off. – whuber Apr 9 '16 at 17:58
• As far as the algorithm goes, there is no difference whatsoever between running it to a specified time or to a specified population, provided you run it long enough. It would still be efficient to run to a much larger than expected population and cut off any values occurring after your target time, rather than running just to your target time, because the algorithm becomes so simple that it parallelizes beautifully (as shown in the R code I posted). – whuber Apr 9 '16 at 18:01
• @whuber: i fixed the simulation to now pass the exponential parameter the inverse of scale (i.e. 1/(size of population * rate)). it still underestimates exponential growth. i understand the parallelization/vectorized version of the code and why it's better but i am not interested in efficiency just in understanding the algorithm, and why my simple version fails. – mvd Apr 9 '16 at 18:26

The whole point of simulation is to show you that such variation is realistic. In fact, there appears to be nothing wrong with your results--except that you haven't yet done enough simulations to appreciate what they are telling you. In this case, the results are especially erratic because the starting population is so small.

Let's run your scenario 500 times up to a population of 300 (rather than a half dozen times to a population of 3):

Just for fun, here's a similar simulation for a population that grows from one individual to its carrying capacity:

I used R for these simulations and plots, because it does one very interesting thing. Since you know in advance that the population will progress in whole steps from the initial population to the final one, you can easily generate the sequence of population values in advance of the simulation. Thus, all that remains is to generate a set of exponentially distributed variates with rates determined by that sequence and accumulate them to simulate the birth times. R performs that with a single command (the line that creates simulation below). It takes about one second. Everything else is just parameter specification and plotting.

(I can get away with such a simple algorithm because I am running these simulations out to a given population target rather than out to a given time endpoint. Obviously the model is the same; all that differs is how I control the length of the simulation.)

rate <- 1
pop.0 <- 1
time.0 <- 0
k <-   0        # Carrying capacity (use 0 or negative when not applicable)
n.final <- 300  # Must not exceed the capacity!
#
# Pre-calculation: populations and the associated rates.
#
n <- pop.0:(n.final-1)
if (k <= 0) r <- rate * n else r <- rate * n * (1 - n/k)
#
# The simulation.
# Each iteration is stored as a column of the result.
#
simulation <- replicate(500, cumsum(c(time.0, rexp(length(n), r))))
#
# Plot the results:
# Set it up, show the overlaid growth curves, then plot a reference curve.
#
plot(range(simulation), c(pop.0, n.final), type="n", ylab="Population", xlab="Time")
apply(simulation, 2, function(x) lines(x, c(n, n.final), col="#00000020"))
if (k <= 0) {curve(pop.0 * exp((x - time.0)*rate), add=TRUE, col="Red", lwd=2)} else
curve(k*(1 - 1/(1+(pop.0/(k-pop.0))*exp(rate*(x-time.0)))), add=TRUE, col="Red", lwd=2)

• It may be worth pointing out this is an unrealistic simulation model for most situations, because it does not allow for the possibility of population reduction. Every one of the simulated curves necessarily is monotonically increasing. The chief reason for using this procedure would seem to be the algorithmic simplification reflected in this one-line R solution. – whuber Apr 8 '16 at 16:42
• thanks for illuminating answer but i don't understand a key part - see my update to post, with new plot, in response to your answer – mvd Apr 9 '16 at 4:42
• @mvd: Rather than adding a substantial amount of new material to your question ---in response to an accepted answer--- it would be better to change it back to what it was, and then ask a separate follow-up question, linking to this question. – Reinstate Monica Apr 18 '18 at 2:41