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))) == 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) expected_answer = math.sqrt(2)- 1 print 'expected probability: %.4f' % (expected_answer) print 'delta from answer: %.4f' % (extinction_probability - expected_answer) if __name__ == '__main__': start = time.clock() main() print 'runtime: %.3f s' % (time.clock() - start) print 'done'