Newton-Pepys, common birthdays and other probability simulations in Python 3
1. Design

I think that using the @mean_of_experiments decorator is not the best approach to this problem:


There might be cases where you want to run the underlying function by itself, for example in order to test it, but the decorator prevents this.
You have to choose the number of experiments, N, when you define the function. This is restrictive: it means that you can't run the function with different values of N.
The only result you get back is the mean. But what if you sometimes want another statistic, such as the standard deviation?


My preferred design would use the principle of separation of concerns. For example, we might write:

def mean_of_experiments(n, f, *args, **kwargs):
    "Return the mean of n calls to f(*args, **kwargs)."
    return mean(f(*args, **kwargs) for _ in range(n))


This meets all my objections:


I can call common_birthday(23) to carry out one experiment, or mean_of_experiments(1000, common_birthday, 23) to find the mean of 1000 experiments.
I can choose a different value of n for each call to mean_of_experiments, I don't have to commit to a particular value in advance.
If I need the standard deviation, I can write std_of_experiments, and this can coexist in the same program with mean_of_experiments.


2. Review


The code is very clear and well documented! In the top 1% of questions on Code Review.
If you're going to use NumPy at all, I think you might as well take advantage of it everywhere, instead of mixing and matching with the Python standard library. For example, in common_birthday you could use numpy.random.randint and numpy.unique:

def common_birthday(k, days=365):
    "Return True iff there are any repeats among k random birthdays."
    return k != np.unique(np.random.randint(1, days, size=k)).size

Similarly, in matching, you could use numpy.random.permutation, numpy.arange, and numpy.any:

def matching(k=52):
    """Return True iff there's any j such that card j appears in the j'th
    place of a shuffled deck of k cards.

    """
    return (np.random.permutation(k) == np.arange(k)).any()

The docstring for boardings doesn't explain the meaning of the scale and N arguments.
The scale argument to boardings seems unnecessary, since the randomly chosen instants will scale along with the arrival times.