Monte Carlo simulation of amoeba population
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 &gt; 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
    num_dead_populations = sum(
        trial(
            amoeba_generations(max_split_per_amoeba),
            lambda n: n == 0,
            lambda n: n &gt;= 10     # Assume a population of &gt;= 10 will never die
        )
        for _ in range(num_trials)
    )
    print("""***starting simulation***
num trials: {num_trials}
max_split_per_amoeba: {max_split_per_amoeba}
extinct outcomes: {num_dead_populations}
extinction probability: {extinction_probability:.4f}
expected probability: {expected_probability:.4f}
delta from answer: {delta:.4f}""".format(
        extinction_probability=num_dead_populations/num_trials,
        expected_probability=2**.5 - 1,
        delta=(num_dead_populations/num_trials) - (2**.5 - 1),
        **locals()
    ))

if __name__ == '__main__':
    main()