# Monte Carlo simulation of amoeba population

I wrote a simple Python simulation to answer the "Amoeba" population question posed here:

A population of amoebas starts with 1. After 1 period that amoeba can divide into 1, 2, 3, or 0 (it can die) with equal probability. What is the probability that the entire population dies out eventually?

The script seems to work fine, but then I started playing with it to investigate curious effects.

If the number of trial generations is too large (num_generations_per_trial), I have problems with performance - the population size gets huge, and the simulation either runs slow or I encounter OverflowError on my brute force FOR loops.

I would appreciate feedback on efficiency options, and also on general code improvements. I know the runs are independent and could be run in parallel. But that is still sort of brute force. I am curious more about making a single thread approach faster.

from __future__ import division
import random
import math
import time

def run_trial(max_split, num_generations):
population = 1
for generation in xrange(num_generations):
for amoeba in xrange(population):
amoeba_split = random.randint(0, max_split)
population -= 1  # remove current amoeba (she will split or die)
population += amoeba_split

if population == 0:
break
return population

def main():
extinct_counter = 0
trials = 10000

max_split_per_amoeba = 3
num_generations_per_trial = 20 # populations can get *massive* as generations increase (memory / overflow errors at 100)

print '***starting simulation***'
print 'num trials: %i' % (trials)
print 'max_split_per_amoeba: %i' % (max_split_per_amoeba)
print 'num_generations_per_trial: %i' % (num_generations_per_trial)

for trial in xrange(trials):
outcome_population = run_trial(max_split_per_amoeba, num_generations_per_trial)
if outcome_population == 0:
extinct_counter += 1

if divmod(trial+1, max(1,int(trials/20)))[1] == 0:
print 'progress: %i trials complete | %i extinction counter | %.4f extinction probability' %  (trial+1, extinct_counter, extinct_counter/(trial+1))

print 'extinct outcomes: %i' % (extinct_counter)
print 'total trials: %i' % (trials)

extinction_probability = extinct_counter / trials
print 'extinction probability: %.4f' % (extinction_probability)

print 'expected probability: %.4f' % (expected_answer)

if __name__ == '__main__':
start = time.clock()
main()
print 'runtime: %.3f s' % (time.clock() - start)
print 'done'

• Overflow errors can be avoided using longs (1L) instead of ints (1); or by switching to python 3 which supports arbitrary large integers. – 301_Moved_Permanently May 5 '16 at 16:19
• @MathiasEttinger - maybe I misunderstood, but I thought the int restriction was built into xrange? "Annoyingly, in Python 2, xrange requires its arguments to fit into a C long" stackoverflow.com/questions/22114088/… – Roberto May 5 '16 at 21:53
• I wasn't aware of that. But anyway, you're using integers, so you hit their limit long before you would have reached it using longs. – 301_Moved_Permanently May 5 '16 at 22:32

## Performance

Limiting each trial to 20 generations is not really helpful. Let's put it this way: the population is expected to increase exponentially with each generation. Once the population reaches, say, 10 amoebas, it's pretty safe to declare that infection has taken hold and the population will never dwindle to 0.

If you run a trial for 20 generations, you can easily get populations in the thousands. The problem with letting populations grow so large is that you have to decide the fate of each individual amoeba at each time step, and that ends up being a lot of work.

If you change num_generations_per_trial to 10, you'll get identical results, and finish much more quickly. But how do you decide what an appropriately low num_generations_per_trial number is? It would be much more intuitive to impose a cap on the population size than on the number of generations.

## Expressiveness

I suggest writing a generator for the population count at each time step. You can easily count elements using the sum() function.

from __future__ import division
from random import randint

def amoeba_generations(max_split=3, population=1):
yield population
while population > 0:
population = sum(randint(0, max_split) for _ in range(population))
yield population

def trial(generator, success_pred, failure_pred):
for g in generator:
if success_pred(g):
return True
if failure_pred(g):
return False
return False

def main():
num_trials = 10000
max_split_per_amoeba = 3
trial(
amoeba_generations(max_split_per_amoeba),
lambda n: n == 0,
lambda n: n >= 10     # Assume a population of >= 10 will never die
)
for _ in range(num_trials)
)
print("""***starting simulation***
num trials: {num_trials}
max_split_per_amoeba: {max_split_per_amoeba}
extinction probability: {extinction_probability:.4f}
expected probability: {expected_probability:.4f}
expected_probability=2**.5 - 1,
**locals()
))

if __name__ == '__main__':
main()

• thank you. I worry that 10 amoeba (or whatever number) might not be valid any more if the other inputs are changed (max_split_per_amoeba, etc). But I understand that 20 generations is equally random, and I see what you are getting at here. It is a good idea and I will certainly think about this a little more. – Roberto May 5 '16 at 22:15
• with performance problems always profiler is huge help. Sometimes your IDE support it, sometimes you must use it directly

• i replace your __main__ section by code

if __name__ == "__main__":
import cProfile # should be at the top of file
cProfile.run('main()')

• when run, i get interesting result

         198906976 function calls in 99.713 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.000    0.000   99.713   99.713 <string>:1(<module>)
1    0.011    0.011   99.713   99.713 ameba.py:20(main)
10000   18.129    0.002   99.697    0.010 ameba.py:7(run_trial)
66292324   59.878    0.000   63.150    0.000 random.py:175(randrange)
66292324   18.417    0.000   81.567    0.000 random.py:238(randint)
10000    0.002    0.000    0.002    0.000 {divmod}
1    0.000    0.000    0.000    0.000 {math.sqrt}
10000    0.002    0.000    0.002    0.000 {max}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
66292324    3.272    0.000    3.272    0.000 {method 'random' of '_random.Random' objects}

• this result means that your program 81.5s of 99s total time is inside random.randint(0, max_split) function. This is important information. You should now only optimize this line, nothing more.

• first though: Can i reduce calls of randint()? Unfortunately, without changing brute-force algorithm it's impossible.

• second though: What does to hell with randrange? I don't use it in my code. The answer: randint calls randrange. So, maybe we call it directly. I replace random.randint(0, max_split) by random.randrange(0, max_split + 1). Result:

         132657778 function calls in 69.706 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.000    0.000   69.706   69.706 <string>:1(<module>)
1    0.010    0.010   69.706   69.706 ameba.py:20(main)
10000   16.422    0.002   69.691    0.007 ameba.py:7(run_trial)
66313887   50.472    0.000   53.269    0.000 random.py:175(randrange)
10000    0.002    0.000    0.002    0.000 {divmod}
1    0.000    0.000    0.000    0.000 {math.sqrt}
10000    0.002    0.000    0.002    0.000 {max}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
66313887    2.797    0.000    2.797    0.000 {method 'random' of '_random.Random' objects}

• it's better, but still this one line takes 53s of 70s total time. We are still need only optimize this line. When i try google this problem, i get similar answers: use numpy. It's not a pure standard library, so it's your choice.

• from my point of view numpy and scipy are "standard libraries" :)

• thanks for this great answer. I want to investigate some more, but I can see already the numpy version does increase speed substantially. (I did note that the random upper bound needs to be increased by 1 for numpy since that max value is not included in numpy's implementation. It is <= in random's version). – Roberto May 5 '16 at 21:51