# Birthday paradox simulation

I'm writing a simulator which checks how many balls we need to occupy every box or how many balls we need for any box has at least 2 balls (birthday paradox). I wrote the script in python and it is so slow. My friend wrote a program in C++ and the script is much faster. For boxes_max_num = 1000 it takes few minutes when in my python script it takes 'eternity'.

There is my code:

import numpy as np
import matplotlib.pyplot as plt

def check_every_box_is_occupied(boxes):
for box in boxes:
if box == 0:
return False
return True

def check_birthday_paradox(boxes):
for box in boxes:
if box >= 2:
return True
return False

def main():
number_of_tests = 250
birthday_paradox_graph = [[], []]
every_box_is_occupied_graph = [[], []]
boxes_max_num = 1000
for number_of_boxes in range(10, boxes_max_num + 1, 1):
print(number_of_boxes)
average_frequency_birthday_paradox = 0
average_frequency_every_box_is_occupied = 0
for index in range(number_of_tests):
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
while True:
boxes[np.random.randint(number_of_boxes)] += 1
if check_birthday_paradox(boxes):
average_frequency_birthday_paradox += number_of_balls
break
number_of_balls += 1
number_of_balls = number_of_boxes
boxes = np.array([0] * number_of_boxes)
while True:
boxes[np.random.randint(number_of_boxes)] += 1
if check_every_box_is_occupied(boxes):
average_frequency_every_box_is_occupied += number_of_balls
break
number_of_balls += 1

plt.rcParams.update({'font.size': 15})
birthday_paradox_graph[0].append(number_of_boxes)
birthday_paradox_graph[1].append(average_frequency_birthday_paradox / number_of_tests)
every_box_is_occupied_graph[0].append(number_of_boxes)
every_box_is_occupied_graph[1].append(average_frequency_every_box_is_occupied / number_of_tests)
plt.figure(1)
plt.plot(birthday_paradox_graph[0], birthday_paradox_graph[1], 'ko')
plt.title("Conajmniej jedna urna ma conajmniej dwie kule")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.figure(2)
plt.title("Wszystkie urny są zapełnione")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.plot(every_box_is_occupied_graph[0], every_box_is_occupied_graph[1], 'ko')
plt.show()

if __name__ == '__main__':
main()


Can you help me to improve the code to make it faster? Is it possible in raw python?

• @Slothario Comments arn't the place to post answers, and con be deleted at any time without notice. – Peilonrayz May 8 '19 at 14:27
• Regarding for number_of_boxes in range(10, boxes_max_num + 1, 1), that's going to perform the following tests for every number between 10 and boxes_max_num +1. Is that the intended behavior? – ShapeOfMatter May 8 '19 at 15:40

## 2 Answers

First, let's break some of the loop bodies out into their own functions. If you're having difficulty reasoning about the behavior of your code, then you can definitely afford the (theoretical) performance hit of a function call.

import numpy as np
import matplotlib.pyplot as plt

def check_every_box_is_occupied(boxes):
for box in boxes:
if box == 0:
return False
return True

def check_birthday_paradox(boxes):
for box in boxes:
if box >= 2:
return True
return False

def run_test(number_of_boxes):
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
result = {
'balls_for_paradox': 0,
'balls_for_full': 0,
}
while True:
boxes[np.random.randint(number_of_boxes)] += 1
if check_birthday_paradox(boxes):
result['balls_for_paradox'] = number_of_balls
break
number_of_balls += 1
number_of_balls = number_of_boxes
boxes = np.array([0] * number_of_boxes)
while True:
boxes[np.random.randint(number_of_boxes)] += 1
if check_every_box_is_occupied(boxes):
result['balls_for_full'] = number_of_balls
break
number_of_balls += 1
return result

def run_tests(number_of_boxes, number_of_tests):
print(number_of_boxes)
average_frequency_birthday_paradox = 0
average_frequency_every_box_is_occupied = 0
for index in range(number_of_tests):
result = run_test(number_of_boxes)
average_frequency_birthday_paradox += result['balls_for_paradox']
average_frequency_every_box_is_occupied += result['balls_for_full']
plt.rcParams.update({'font.size': 15})
return {
'average_frequency_birthday_paradox': average_frequency_birthday_paradox / number_of_tests,
'average_frequency_every_box_is_occupied': average_frequency_every_box_is_occupied / number_of_tests,
}

def main():
number_of_tests = 250
birthday_paradox_graph = [[], []]
every_box_is_occupied_graph = [[], []]
boxes_max_num = 1000
for number_of_boxes in range(10, boxes_max_num + 1, 1):
results = run_tests(number_of_boxes, number_of_tests)
birthday_paradox_graph[0].append(number_of_boxes)
birthday_paradox_graph[1].append(results['average_frequency_birthday_paradox'])
every_box_is_occupied_graph[0].append(number_of_boxes)
every_box_is_occupied_graph[1].append(results['average_frequency_every_box_is_occupied'])
plt.figure(1)
plt.plot(birthday_paradox_graph[0], birthday_paradox_graph[1], 'ko')
plt.title("Conajmniej jedna urna ma conajmniej dwie kule")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.figure(2)
plt.title("Wszystkie urny są zapełnione")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.plot(every_box_is_occupied_graph[0], every_box_is_occupied_graph[1], 'ko')
plt.show()

if __name__ == '__main__':
main()


The above is a pretty literal translation of your code. I introduced some dictionaries as return values; tuples would have worked but would have required some kind of documentation. Explicit classes would be more maintainable, but let's not worry about that for this project.

Next, I'm pretty sure you have a math error in your calculations of the "every box occupied" condition. The line where you reset number of balls should be resetting it to 1, right?
I'm also going to use explicit else statements in this round, because I think they look better.

def run_test(number_of_boxes):
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
result = {
'balls_for_paradox': 0,
'balls_for_full': 0,
}
while True:
boxes[np.random.randint(number_of_boxes)] += 1
if check_birthday_paradox(boxes):
result['balls_for_paradox'] = number_of_balls
break
else:
number_of_balls += 1
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
while True:
boxes[np.random.randint(number_of_boxes)] += 1
if check_every_box_is_occupied(boxes):
result['balls_for_full'] = number_of_balls
break
else:
number_of_balls += 1
return result


Now we can implement Slothario's suggestion, but we'll have to rethink how we're running these tests. The birthday paradox is easy enough, but to avoid checking every cell for the "all occupied" condition, we need to remember cells we've already visited. We can think of this as crossing items off a list.

def run_test(number_of_boxes):
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
result = {
'balls_for_paradox': 0,
'balls_for_full': 0,
}
while True:
box = np.random.randint(number_of_boxes)
boxes[box] += 1
if 2 <= boxes[box]:
result['balls_for_paradox'] = number_of_balls
break
else:
number_of_balls += 1
number_of_balls = 1
boxes = set(range(number_of_boxes))
while True:
box = np.random.randint(number_of_boxes)
if box in boxes:
boxes.remove(box)
if not boxes:
result['balls_for_full'] = number_of_balls
break
else:
number_of_balls += 1
return result


We have two similar loops. In a lot of contexts I would suggest that they're different enough that they should go in completely separate functions, but since we're going for speed, maybe we should combine them.

def run_test(number_of_boxes):
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
unoccupied_indexes = set(range(number_of_boxes))
result = {
'balls_for_paradox': 0,
'balls_for_full': 0,
}
while not (result['balls_for_paradox'] and result['balls_for_full']):
box = np.random.randint(number_of_boxes)
if not result['balls_for_paradox']:
boxes[box] += 1
if 2 <= boxes[box]:
result['balls_for_paradox'] = number_of_balls
if not result['balls_for_full']:
if box in unoccupied_indexes:
unoccupied_indexes.remove(box)
if not unoccupied_indexes:
result['balls_for_full'] = number_of_balls
number_of_balls += 1
return result


Let's tidy up some other little things. You had a call to set the font-size of your plot inside one of the loops, I don't think that belonged. I also removed a print statement, which takes up more time than you might think. I'm also replacing several of the loops with list comprehensions, which won't improve performance but it's a little nicer to read. We could go a lot further in this direction, but I wanted to keep the basic structure of your code intact.

import numpy as np
import matplotlib.pyplot as plt

def check_every_box_is_occupied(boxes):
for box in boxes:
if box == 0:
return False
return True

def check_birthday_paradox(boxes):
for box in boxes:
if box >= 2:
return True
return False

def run_test(number_of_boxes):
number_of_balls = 1
boxes = np.array([0] * number_of_boxes)
unoccupied_indexes = set(range(number_of_boxes))
result = {
'balls_for_paradox': 0,
'balls_for_full': 0,
}
while not (result['balls_for_paradox'] and result['balls_for_full']):
box = np.random.randint(number_of_boxes)
if not result['balls_for_paradox']:
boxes[box] += 1
if 2 <= boxes[box]:
result['balls_for_paradox'] = number_of_balls
if not result['balls_for_full']:
if box in unoccupied_indexes:
unoccupied_indexes.remove(box)
if not unoccupied_indexes:
result['balls_for_full'] = number_of_balls
number_of_balls += 1
return result

def run_tests(number_of_boxes, number_of_tests):
results = [run_test(number_of_boxes) for _ in range(number_of_tests)]
return {
'average_frequency_birthday_paradox': sum([r['balls_for_paradox'] for r in results]) / number_of_tests,
'average_frequency_every_box_is_occupied': sum([r['balls_for_full'] for r in results]) / number_of_tests,
}

def main():
number_of_tests = 250
boxes_max_num = 1000
all_results = [{
'n': number_of_boxes,
'results': run_tests(number_of_boxes, number_of_tests)
}
for number_of_boxes
in range(10, boxes_max_num + 1, 1)]
birthday_paradox_graph = [
[r['n'] for r in all_results],
[r['results']['average_frequency_birthday_paradox'] for r in all_results]
]
every_box_is_occupied_graph = [
[r['n'] for r in all_results],
[r['results']['average_frequency_every_box_is_occupied'] for r in all_results]
]
plt.rcParams.update({'font.size': 15})
plt.figure(1)
plt.plot(birthday_paradox_graph[0], birthday_paradox_graph[1], 'ko')
plt.title("Conajmniej jedna urna ma conajmniej dwie kule")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.figure(2)
plt.title("Wszystkie urny są zapełnione")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.plot(every_box_is_occupied_graph[0], every_box_is_occupied_graph[1], 'ko')
plt.show()

if __name__ == '__main__':
main()


*Runs in ...
my paitence ran out. But if I follow Pelionrayz and use boxes_max_num=1000 and number_of_tests=1, then it takes about 7 seconds.

(just a heads up that I didn't test ever iteration of this; so i don't know if the earlier versions will actually run.)

• An embarrassing error on the "code review" SE: I left in the two unused test functions! – ShapeOfMatter May 8 '19 at 17:55

# Format

Firstly when performance is a concern you should firstly make the code as readable as possible. This is as when looking for performance gains it reduces the readability of the code, and as this is an incremental process it slowly degrades your code to possibly unmanageable levels, whilst making it harder for you to accurately improve performance. And is basically a modern explanation of the house of the rock parable.

You have some issues a linter can notice:

• Imports should be ordered alphabetically so they are easier to read when there's lots of them.
• index is unused, and so should be named _ or _index.
• You should try to keep lines no longer than 79 characters wide.
• Don't mix quote delimiters, pick either ' or ".
• It's good that you have kept PEP8 spacing between your functions, it should be noted this spacing is between functions and anything, and so having one space between main and the main guard is a violation of this.

From here I'd also:

• Split main into a program that performs the birthday_problem and the one that is main.
• Replace check_every_box_is_occupied with all.
• Replace check_birthday_paradox with any and a comprehension.
• Move the generation of infinite amounts of random numbers out of birthday_problem, by creating an infinite generator.
• Wrap the infinite generator in an enumerate to simplify the number_of_balls incrementation.
• Move the while True loops into their own functions.
• Change the for index loop into two comprehensions wrapped in a sum.
• Simplify your naming, you have very long variables that only make the code harder to understand.
• It'd be simpler if you set birthday_paradox_graph[0] to the same range that you iterate over.

And so would get the following code. I have left main uncleaned:

import matplotlib.pyplot as plt
import numpy as np

def random_numbers(limit):
while True:
yield np.random.randint(limit)

def check_birthday_paradox(values):
return any(value >= 2 for value in values)

def simulate_one_pairing(size):
boxes = np.array([0] * size)
for number_of_balls, choice in enumerate(random_numbers(size), 1):
boxes[choice] += 1
if check_birthday_paradox(boxes):
return number_of_balls

def simulate_all_days_set(size):
boxes = np.array([0] * size)
for number_of_balls, choice in enumerate(random_numbers(size), size):
boxes[choice] += 1
if all(boxes):
return number_of_balls

def birthday_problem(tests, boxes_limit):
domain = range(10, boxes_limit + 1, 1)
paired = [domain, []]
every_box_is_occupied_graph = [domain, []]
for boxes in domain:
total = sum(simulate_one_pairing(boxes) for _ in range(tests))
paired[1].append(total / tests)

total = sum(simulate_all_days_set(boxes) for _ in range(tests))
every_box_is_occupied_graph[1].append(total / tests)
return paired, every_box_is_occupied_graph

def main():
number_of_tests = 250
boxes_max_num = 1000
birthday_paradox_graph, every_box_is_occupied_graph = birthday_problem(number_of_tests, boxes_max_num)
plt.rcParams.update({'font.size': 15})
plt.figure(1)
plt.plot(birthday_paradox_graph[0], birthday_paradox_graph[1], 'ko')
plt.title("Conajmniej jedna urna ma conajmniej dwie kule")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.figure(2)
plt.title("Wszystkie urny są zapełnione")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.plot(
every_box_is_occupied_graph[0],
every_box_is_occupied_graph[1],
'ko')
plt.show()

if __name__ == '__main__':
main()


From this I can see a couple of problems off the bat.

1. Why does simulate_all_days_set's start enumerate start at size?
I'll assume this is a mistake.
2. The function check_birthday_paradox is wasting time as it runs in $$\O(n)\$$ time, where all you need to check is that boxes[choice] is equal or greater than two.
3. You're building two simulations each test simulate_one_pairing and simulate_all_days_set, you can just make this a single loop and reduce duplicate simulations.
4. The performance of np.random.randint should be tested against when passing size=None and size=n. To see if chunking can reduce the performance.

And yes it does increase performance: <source>

And so I'd use chunk at the same size as size. You may be able to get better performance if you use a different value however.

import matplotlib.pyplot as plt
import numpy as np
import timeit

def random_numbers(limit):
while True:
yield from np.random.randint(limit, size=limit)

def simulate(size):
boxes = np.array([0] * size)
pair = None
all_set = None
for iteration, choice in enumerate(random_numbers(size), 1):
boxes[choice] += 1
if pair is None and boxes[choice] >= 2:
pair = iteration
if all_set is not None:
break
if all_set is None and all(boxes):
all_set = iteration
if pair is not None:
break
return pair, all_set

def birthday_problem(tests, boxes_limit):
domain = range(10, boxes_limit + 1, 1)
paired = [domain, []]
all_set = [domain, []]
for boxes in domain:
pairs, all_sets = zip(*(simulate(boxes) for _ in range(tests)))
paired[1].append(sum(pairs) / tests)
all_set[1].append(sum(all_sets) / tests)
return paired, all_set

def main():
start = timeit.default_timer()
number_of_tests = 20
boxes_max_num = 200
birthday_paradox_graph, every_box_is_occupied_graph = birthday_problem(number_of_tests, boxes_max_num)
print(timeit.default_timer() - start)
plt.rcParams.update({'font.size': 15})
plt.figure(1)
plt.plot(birthday_paradox_graph[0], birthday_paradox_graph[1], 'ko')
plt.title("Conajmniej jedna urna ma conajmniej dwie kule")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.figure(2)
plt.title("Wszystkie urny są zapełnione")
plt.xlabel("Liczba urn")
plt.ylabel("Średnia liczba kul potrzebna do spełnienia warunku")
plt.plot(
every_box_is_occupied_graph[0],
every_box_is_occupied_graph[1],
'ko')
plt.show()

if __name__ == '__main__':
main()


I tested both yours and mine and at number_of_tests=20 and boxes_max_num=200. You can see this in the code above. Yours unchanged runs in ~13.5-14 seconds. Mine runs in ~3.5-4 seconds.

My code takes ~16.5 seconds to run at boxes_max_num=1000 and number_of_tests=1. This means it will take roughly an hour to run all of your tests. As I've taken the low-hanging fruit if you need it to be faster, then you'll have to use a profiler. Your code also is very unlikely to be as fast as your friends, as Python is slow.