Simulating the Birthday Problem

The code that does all the simulation and plotting is listed below.

import numpy as np
import matplotlib.pylab as plt

def calculate(room = 20):
    return 1 - np.prod((365 - np.arange(room))/365)

def simulate(room = 20, n_sim = 1000):
    hit = 0
    for i in range(n_sim):
        n_bday = np.unique(np.random.randint(1, 365, room)).shape[0]
        hit += n_bday == room
    return 1 - hit/n_sim

plt.plot([calculate(r) for r in range(1, 35)], label="calculated")
plt.plot([simulate(r) for r in range(1, 35)], label="simulated")
plt.legend();