# Simulation of an alien population

### 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?

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)


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.

• Depending on the platform (and especially in python) working with native machine words is faster even if it seems wasteful to only use 2-3 bits out of 64 or 32. – Aaron Dec 10 '18 at 15:58