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

Question

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) >= 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:

1. 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%.
2. 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
3. dice: Problem posed by a gambler to Newton; dice(6, 1) > dice(12, 2) > dice(18, 3)
4. 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?

Answer 1

functools.wraps

Unrelated to your main question, I advice you to use functools.wraps. This way, your methods metadata get transferred to the wrapped method. For example dice.__doc__ returns None without it, but the doctstring if you include it.

repetitions

It is possible to programatically change the number of repetitions per function. If you look at functools.lru_cache, where you need to specify the size of the cache too, you can see you just need an extra level of wrapper function:

def mean_of_experiments(N=100_000):
    def inner_decorator(experiment_func):
        '''Decorator to repeat any Bernoulli trial N times and return probability of success'''
        @wraps(experiment_func)
        def wrapper(*args, **kwargs):
            # print(f"{N} repititions")
            return round(mean(experiment_func(*args, **kwargs) for _ in range(N)), 3)
        return wrapper
    return inner_decorator

where the inner_decorator is the old mean_of_experiments method

Dynamically set repetitions

A decorator is nothing but syntactic sugar for decorator(func), so you can do it by not decorating the functions themselves, but doing the 'decorating' when calling them:

mean_of_experiments(N=100)(dice(100, 20))

Dynamically set repetitions 2:

Another approach is to pass the repetitions on with the kwargs to the experiment_func

def mean_of_experiments_2(experiment_func):
    
    '''Decorator to repeat any Bernoulli trial N times and return probability of success'''
    def wrapper(*args, **kwargs):
        repetitions = kwargs.pop('repetitions', 100_000)
        # print(f"{repetitions} repetitions")
        return round(mean(experiment_func(*args, **kwargs) for _ in range(repetitions)), 3)
    return wrapper

@mean_of_experiments_2
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) >= c

and then call it like this: dice(6, 4, repetitions=100)

The main caveat here is not to pick repetitions as argument to any of the experiment_funcs

Further remarks

_ in ints

you wrote 100000 in 2 ways: 100000 and 1_00_000. The second way looks incorrect to me (it is syntactically correct, but I would group numbers per 3). The way I would write it is 100_000. You can pick another way, but stay consistent

common_birthday

since the rands as you make it will always contain k elements, it would be more efficient to skip this list, and immediately make the set:

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) != k

dice

For large n, this can become a very long list. An alternative approach would be to either use a collections.Counter, or use sum

return  sum(randint(1, 6) == 6 for _ in range(n)) >= c

or

return Counter(randint(1, 6) for _ in range(n))[6] >= c

boardings

If, instead of depending on numpy you want a native implementations, you can use itertools.accumulate and random.expovariate.

def boardings_native(scale=5.0, N=100_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 = list(accumulate(expovariate(lambd=1/scale) for _ in range(N)))

    @mean_of_experiments(3)
    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()

If I look at the performance, it is even faster than the numpy implementation, but that is probably because generating the arrivals is an iteration instead of a vectorised operation

%timeit boardings_native()

52 ms ± 1.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit boardings()

143 ms ± 713 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Answer 2

1. Design

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

1. 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.

2. 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.

3. 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:

1. 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.

2. 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.

3. 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

1. The code is very clear and well documented! In the top 1% of questions on Code Review.

2. 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
   

3. 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()
   

4. The docstring for boardings doesn't explain the meaning of the scale and N arguments.

5. The scale argument to boardings seems unnecessary, since the randomly chosen instants will scale along with the arrival times.

Answer 3

There is a function choices in random library that returns a list of elements chosen from the population with replacement. It could simplify the methods common_birthday and dice. Though, I believe, this decreases readability by a little.

from random import choices

def common_birthday(k):
    return len(set(choices(range(1, 366), k=k))) != k

def dice(n, c):
    return choices(range(1, 7), k=n).count(6) >= c

I removed the docstring and decorator for exposition of business logic.