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 &gt;= 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?
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 &gt;= 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 &lt;= 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 &lt;= 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 &gt;= 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 &lt;= 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.)