Simulation of an alien population

Question

Background

I've come across a puzzle (Problem 10 from this list of sample interview questions):

One day, an alien comes to Earth. Every day, each alien does one of four things, each with equal probability to:

• Kill himself
• Do nothing
• Split himself into two aliens (while killing himself)
• split himself into three aliens (while killing himself)

What is the probability that the alien species eventually dies out entirely?

Unfortunately, I haven't been able to solve the problem theoretically. Then I moved on to simulate it with a basic Markov Chain and Monte Carlo simulation in mind.

This was not asked to me in an interview. I learned the problem from a friend, then found the link above while searching for mathematical solutions.

Reinterpreting the question

We start with the number of aliens n = 1. n has a chance to not change, be decremented by 1, be incremented by 1, and be incremented by 2, %25 for each. If n is incremented, i.e. aliens multiplied, we repeat this procedure for n times again. This corresponds to each alien will do its thing again. I have to put an upper limit though, so that we stop simulating and avoid a crash. n is likely to increase and we're looping n times again and again.

If aliens somehow go extinct, we stop simulating again since there's nothing left to simulate.

After n reaches zero or the upper limit, we also record the population (it will be either zero or some number >= max_pop).

I repeat this many times and record every result. At the end, number of zeros divided by total number of results should give me an approximation.

The code

from random import randint
import numpy as np

pop_max = 100
iter_max = 100000

results = np.zeros(iter_max, dtype=int)

for i in range(iter_max):
    n = 1
    while n > 0 and n < pop_max:
        for j in range(n):
            x = randint(1, 4)
            if x == 1:
                n = n - 1
            elif x == 2:
                continue
            elif x == 3:
                n = n + 1
            elif x == 4:
                n = n + 2
    results[i] = n

print( np.bincount(results)[0] / iter_max )

iter_max and pop_max can be changed indeed, but I thought if there is 100 aliens, the probability that they go extinct would be negligibly low. This is just a guess though, I haven't done anything to calculate a (more) proper upper limit for population.

This code gives promising results, fairly close to the real answer which is approximately %41.4.

Some outputs

> python aliens.py
0.41393
> python aliens.py
0.41808
> python aliens.py
0.41574
> python aliens.py
0.4149
> python aliens.py
0.41505
> python aliens.py
0.41277
> python aliens.py
0.41428
> python aliens.py
0.41407
> python aliens.py
0.41676

Aftermath

I'm okay with the results but I can't say the same for the time this code takes. It takes around 16-17 seconds :)

How can I improve the speed? How can I optimize loops (especially the while loop)? Maybe there's a much better approach or better models?

Answer 1

There are some obvious possible simplifications for elegance, if not necessarily for speed.

The while condition should be written as a double-ended inequality:

while 0 < n < pop_max:
    …

The variable j is unused. The convention is to use _ as the name of a "throwaway" variable.

The if-elif chain can be eliminated with a smarter randint() call:

for j in range(n):
    n += randint(-1, 2)

NumPy is overkill here, when all you want to know whether the population went extinct. The built-in sum() function can do the counting for you.

Each simulation run is independent. I'd put the code in a function for readability.

from random import randint

def population(pop_max=100):
    n = 1
    while 0 < n < pop_max:
        for _ in range(n):
            n += randint(-1, 2)
    return n

iterations = 100000
print(sum(population() == 0 for _ in range(iterations)) / iterations)

Answer 2

numpy has a function to generate a random array, this might be faster than generating a random number within the inner loop. https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.randint.html

You can also try generating a larger 32 or 64-bit number, and shifting and masking the whole time to get 2 random bits. However, this is a bit far-fetched.