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:
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/edice
: 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 ofscale
.
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?