Newton-Pepys, common birthdays and other probability simulations in Python 3 I have been learning Probability from Statistics 110 and trying to simulate problems that have an unintuitive answer

from numpy import cumsum
from statistics import mean
from numpy.random import exponential
from random import randint, sample, uniform
from bisect import bisect_left

def mean_of_experiments(experiment_func, N=100000):
    '''Decorator to repeat any Bernoulli trial N times and return probability of success'''
    def wrapper(*args, **kwargs):
        return round(mean(experiment_func(*args, **kwargs) for _ in range(N)), 3)
    return wrapper

@mean_of_experiments
def common_birthday(k):
    '''Simulates an experiment to generate k independent uniformly random birthdays and check if there are any repeat birthdays'''
    rands = [randint(1, 365) for _ in range(k)]
    return len(rands) != len(set(rands))

@mean_of_experiments
def matching(k=52):
    '''Simulates an experiment to permute 'k' cards and check if any jth card's value is j'''
    idx_labels = enumerate(sample(range(k), k))
    return any(idx == label for idx, label in idx_labels)

@mean_of_experiments
def dice(n, c):
    '''Simulates an experiment to roll 'n' dice and and check if count of 6's is at least c'''
    return [randint(1, 6) for _ in range(n)].count(6) &gt;= c

def boardings(scale=5.0, N=1_00_000):
    '''Simulates an experiment where arrival of buses at stop follows a Poisson process and finds avg. inter-arrival time at a random instant'''
    arrivals = cumsum([exponential(scale=scale) for _ in range(N)])

    @mean_of_experiments
    def wait():
        boarding_idx = bisect_left(arrivals, uniform(0, arrivals[-1]))
        missed_bus = 0 if boarding_idx == 0 else arrivals[boarding_idx - 1]
        return arrivals[boarding_idx] - missed_bus

    return wait()


Probability Problems:


common_birthday: Given k people, what's the probability that any 2 people share a birthday? For 23 people, it's >50% and for 50 people, >97%.
matching: Given 52 cards with a unique label in 0...n-1, they are shuffled. What's the probability that any j^th card's label is j? Answer's very close to 1-1/e
dice: Problem posed by a gambler to Newton; dice(6, 1) > dice(12, 2) > dice(18, 3)
boardings: Given that a bus company's buses follow a Poisson process, if a person visits a stop at some random time, what was that inter-arrival time? From Tsitsiklis' lectures; answer is 2*scale for any value of scale.


What I wish I could do in my code is dynamically set the value of N in the decorator mean_of_experiments. Is that possible?
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.