#!/usr/bin/env python
# coding: utf-8
# # Table of Contents
#
# # Bernoulli and binomial distribution
# - References: [Bernoulli distribution on Wikipedia](https://en.wikipedia.org/wiki/Bernoulli_distribution) and [Binomial distribution on Wikipedia](https://en.wikipedia.org/wiki/Binomial_distribution#Generating_binomial_random_variates).
#
# The Bernoulli distribution of mean $p\in[0,1]$ is defined as the distribution on $\{0,1\}$ such that $\mathbb{P}(X=1) = p$ and $\mathbb{P}(X=0) = 1-p$.
#
# If $X$ follows a Binomial distribution of mean $p\in[0,1]$ and $n$ samples, $X$ is defined as the sum of $n$ independent and identically distributed (iid) samples from a Bernoulli distribution of mean $p$, that is $X\in\{0,\dots,n\}$ ($X\in\mathbb{N}$) and $\forall k\in\{0,\dots,n\}, \mathbb{P}(X=k) = {n \choose k} p^k (1-p)^{n-k}$.
# ## Requirements
# Let's import the modules required for this notebook.
# In[1]:
import numpy as np
import matplotlib.pyplot as plt
# In[2]:
get_ipython().run_line_magic('load_ext', 'cython')
# In[3]:
get_ipython().run_line_magic('load_ext', 'watermark')
get_ipython().run_line_magic('watermark', '-a "Lilian Besson (Naereen)" -i -v -p numpy,matplotlib,cython')
# ## A naive generator
# Using the pseudo-random generator of (`float`) random numbers in $[0,1]$ from the `random` or `numpy.random` module, we can easily generate a sample from a Bernoulli distribution.
# In[5]:
import random
def uniform_01() -> float:
return random.random()
# In[7]:
[ uniform_01() for _ in range(5) ]
# It's very quick now:
# In[12]:
def bernoulli(p: float) -> int:
return 1 if uniform_01() <= p else 0
# In[13]:
[ bernoulli(0) for _ in range(5) ]
# In[14]:
[ bernoulli(0.12345) for _ in range(5) ]
# In[15]:
[ bernoulli(1) for _ in range(5) ]
# So we can naively generate samples from a Binomial distribution by summing iid samples generated using this `bernoulli` function.
# In[16]:
def naive_binomial(n: int, p: float) -> int:
result = 0
for k in range(n): # sum of n iid samples from Bernoulli(p)
result += bernoulli(p)
return result
# For example :
# In[17]:
[ naive_binomial(10, 0.1) for _ in range(5) ]
# In[18]:
[ naive_binomial(10, 0.5) for _ in range(5) ]
# In[19]:
[ naive_binomial(10, 0.9) for _ in range(5) ]
# We can quickly illustrate the generated distribution, to check it has the correct "shape":
# In[21]:
m = 1000
n = 10
p = 0.12345
X = [ naive_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[22]:
m = 1000
n = 10
p = 0.5
X = [ naive_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[23]:
m = 1000
n = 10
p = 0.98765
X = [ naive_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# ## The generator included in `numpy.random`
# In[28]:
def numpy_binomial(n: int, p: float) -> int:
return np.random.binomial(n, p)
# Let's try this out:
# In[29]:
[ numpy_binomial(10, 0.1) for _ in range(5) ]
# In[30]:
[ numpy_binomial(10, 0.5) for _ in range(5) ]
# In[31]:
[ numpy_binomial(10, 0.9) for _ in range(5) ]
# Let's plot this out also.
# In[32]:
m = 1000
n = 10
p = 0.12345
X = [ numpy_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[33]:
m = 1000
n = 10
p = 0.5
X = [ naive_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[34]:
m = 1000
n = 10
p = 0.98765
X = [ naive_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# ## An efficient generator using the inverse transform method
#
# 1. We start by computing the binomial coefficients and then the probability $\mathbb{P}(X=k)$ for $k\in\{0,\dots,n\}$, if $X\sim Bin(n, p)$, and,
# 2. Then use this to write a generator of Binomial-distributed random values.
# 3. This function is then simplified to inline all computations.
# 4. We propose a fast and simple Cython implementation, to be as efficient as possible, and hopefully comparably efficient when compared against the implementation in Numpy.
# ### Explicit computation of the probabilities
# In[36]:
def binomial_coefficient(n: int, k: int) -> int:
"""From https://en.wikipedia.org/wiki/Binomial_coefficient#Binomial_coefficient_in_programming_languages"""
if k < 0 or k > n:
return 0
if k == 0 or k == n:
return 1
k = min(k, n - k) # take advantage of symmetry
c = 1
for i in range(k):
c = (c * (n - i)) / (i + 1)
return c
# In[37]:
def proba_binomial(n: int, p: float, k: int) -> float:
"""Compute {n \choose k} p^k (1-p)^(n-k)"""
q = 1.0 - p
return binomial_coefficient(n, k) * p**k * q**(n-k)
# ### First function using the inversion method
# This first function is a generic implementation of the discrete inverse transform method.
# For more details, see [the Wikipedia page](https://en.wikipedia.org/wiki/Inverse_transform_sampling).
#
# > Inverse transformation sampling takes uniform samples of a number $u$ between $0$ and $1$, interpreted as a probability, and then returns the largest number $x$ from the domain of the distribution $\mathbb{P}(X)$ such that $\mathbb{P}(-\infty int:
probas = [ compute_proba(x) for x in range(xmin, xmax + 1) ]
result = xmin
current_proba = 0
one_uniform_sample = uniform_01()
while current_proba <= one_uniform_sample:
current_proba += probas[result]
result += 1
return result - 1
# In[48]:
def first_inversion_binomial(n: int, p: float) -> int:
def compute_proba(x):
return proba_binomial(n, p, x)
xmax = n
xmin = 0
return inversion_method(compute_proba, xmax, xmin=xmin)
# Let's try out.
# In[51]:
[ first_inversion_binomial(10, 0.1) for _ in range(5) ]
# In[52]:
[ first_inversion_binomial(10, 0.5) for _ in range(5) ]
# In[53]:
[ first_inversion_binomial(10, 0.9) for _ in range(5) ]
# It seems to work as wanted!
# ### Simplified code of the inversion method
#
# The previous function as a few weaknesses: it stores the $n+1$ values of $\mathbb{P}(X=k)$ before hand, it computes all of them even if the `for` loop of the inversion method stops in average before the end (in average, it takes $np$ steps, which can be much smaller than $n$ for small $p$).
# Furthermore, the computations of both the binomial coefficients and the values $p^k (1-p)^{n-k}$ is using powers and not iterative multiplications, leading to more rounding errors.
#
# We can solve all these issues by inlining all the computations.
# In[104]:
def inversion_binomial(n: int, p: float) -> int:
if p <= 1e-10:
return 0
if p >= 1 - 1e-10:
return n
if p > 0.5: # speed up by computing for q and then substracting
return n - inversion_binomial(n, 1.0 - p)
result = 0
q = 1.0 - p
current_proba = q**n
cum_proba = current_proba
one_uniform_sample = uniform_01()
while cum_proba <= one_uniform_sample:
current_proba *= (p * (n - result)) / (q * (result + 1))
cum_proba += current_proba
result += 1
return result
# Let's try out.
# In[64]:
[ inversion_binomial(10, 0.1) for _ in range(5) ]
# In[82]:
[ inversion_binomial(10, 0.5) for _ in range(5) ]
# In[67]:
[ inversion_binomial(10, 0.9) for _ in range(5) ]
# It seems to work as wanted!
#
# And now the storage is indeed $O(1)$, and the computation time is $O(x)$ if the return value is $x$, so the mean computation time is $O(np)$.
#
# Note that if $p=1/2$, then $O(np) = O(n/2) = O(n)$, and thus this improved method using the inversion method is (asymptotically) as costly as the naive method (the first method which consists of summing $n$ iid samples from a Bernoulli of mean $p$).
# Let's plot this out also.
# In[68]:
m = 1000
n = 10
p = 0.12345
X = [ inversion_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[87]:
m = 1000
n = 10
p = 0.5
X = [ inversion_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[88]:
m = 1000
n = 10
p = 0.98765
X = [ inversion_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# ### In Cython
# In[89]:
get_ipython().run_line_magic('load_ext', 'cython')
# In[112]:
get_ipython().run_cell_magic('cython', '--annotate', '\nimport random\n\ndef cython_inversion_binomial(int n, double p) -> int:\n if p <= 1e-9:\n return 0\n if p >= 1 - 1e-9:\n return n\n if p > 0.5: # speed up by computing for q and then substracting\n return n - cython_inversion_binomial(n, 1.0 - p)\n cdef int result = 0\n cdef double q = 1.0 - p\n cdef double current_proba = q**n\n cdef double cum_proba = current_proba\n cdef double one_uniform_sample = random.random()\n while cum_proba < one_uniform_sample:\n current_proba *= (p * (n - result)) / (q * (result + 1))\n cum_proba += current_proba\n result += 1\n return result\n')
# Let's try out.
# In[113]:
[ cython_inversion_binomial(10, 0.1) for _ in range(5) ]
# In[114]:
[ cython_inversion_binomial(10, 0.5) for _ in range(5) ]
# In[115]:
[ cython_inversion_binomial(10, 0.9) for _ in range(5) ]
# It seems to work as wanted!
# Let's plot this out also.
# In[116]:
m = 1000
n = 10
p = 0.12345
X = [ cython_inversion_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[117]:
m = 1000
n = 10
p = 0.5
X = [ cython_inversion_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# In[118]:
inversion_binomialm = 1000
n = 10
p = 0.98765
X = [ cython_inversion_binomial(n, p) for _ in range(m) ]
plt.figure()
plt.hist(X)
plt.title(f"{m} samples from a Binomial distribution with n = {n} and p = {p}.")
plt.show()
# ## Numerical experiments to check time cost of the different versions
# In[119]:
n = 100
# In[120]:
naive_binomial
first_inversion_binomial
inversion_binomial
cython_inversion_binomial
numpy_binomial
# We can use the `%timeit` magic to check the (mean) computation time of all the previously mentioned functions:
# In[106]:
get_ipython().run_line_magic('timeit', 'naive_binomial(n, 0.123456)')
get_ipython().run_line_magic('timeit', 'first_inversion_binomial(n, 0.123456)')
get_ipython().run_line_magic('timeit', 'inversion_binomial(n, 0.123456)')
get_ipython().run_line_magic('timeit', 'cython_inversion_binomial(n, 0.123456)')
get_ipython().run_line_magic('timeit', 'numpy_binomial(n, 0.123456)')
# Apparently, our `cython` method is faster than the function from `numpy`!
#
# We also check that our first naive implementation of the inversion method was suboptimal, as announced, because of its pre computation of all the values of $\mathbb{P}(X=k)$.
# However, we check that the naive method, using the sum of $n$ binomial samples, is as comparably efficient to the pure-Python inversion-based method (for this small $n=100$).
# In[121]:
get_ipython().run_line_magic('timeit', 'naive_binomial(n, 0.5)')
get_ipython().run_line_magic('timeit', 'first_inversion_binomial(n, 0.5)')
get_ipython().run_line_magic('timeit', 'inversion_binomial(n, 0.5)')
get_ipython().run_line_magic('timeit', 'cython_inversion_binomial(n, 0.5)')
get_ipython().run_line_magic('timeit', 'numpy_binomial(n, 0.5)')
# In[108]:
get_ipython().run_line_magic('timeit', 'naive_binomial(n, 0.987654)')
get_ipython().run_line_magic('timeit', 'first_inversion_binomial(n, 0.987654)')
get_ipython().run_line_magic('timeit', 'inversion_binomial(n, 0.987654)')
get_ipython().run_line_magic('timeit', 'cython_inversion_binomial(n, 0.987654)')
get_ipython().run_line_magic('timeit', 'numpy_binomial(n, 0.987654)')
# It's quite awesome to see that our inversion-based method is more efficient that the numpy function, both in the pure-Python and the Cython versions!
# But it's weird, as the numpy function is... based on the inversion method, and itself written in C!
#
# > See the source code, [numpy/distributions.c line 426](https://github.com/numpy/numpy/blob/7c41164f5340dc998ea1c04d2061f7d246894955/numpy/random/mtrand/distributions.c#L426) (on the 28th February 2019, commit 7c41164).
#
# But the trick is that the implementation in numpy uses the inversion method (running in $\Omega(np)$) if $pn < 30$, and a method denoted "BTPE" otherwise.
# I need to work on this method! The BTPE algorithm is much more complicated, and it is described in the following paper:
#
# > Kachitvichyanukul, V.; Schmeiser, B. W. (1988). "Binomial random variate generation". Communications of the ACM. 31 (2): 216–222. [doi:10.1145/42372.42381](https://doi.org/10.1145%2F42372.42381).
#
# > See the source code, [numpy/distributions.c line 263](https://github.com/numpy/numpy/blob/7c41164f5340dc998ea1c04d2061f7d246894955/numpy/random/mtrand/distributions.c#L263) (on the 28th February 2019, commit 7c41164).
# ## Checking that sampling from $Bin(n,p)$ requires a time $\Omega(n)$.
# In[123]:
n = 100
get_ipython().run_line_magic('timeit', 'naive_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'inversion_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'cython_inversion_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'numpy_binomial(n, random.random())')
# In[124]:
n = 1000
get_ipython().run_line_magic('timeit', 'naive_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'inversion_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'cython_inversion_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'numpy_binomial(n, random.random())')
# In[ ]:
n = 10000
get_ipython().run_line_magic('timeit', 'naive_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'inversion_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'cython_inversion_binomial(n, random.random())')
get_ipython().run_line_magic('timeit', 'numpy_binomial(n, random.random())')
# As we can see, our inversion method (no matter the implementation) runs in $O(n)$ (for $p$ in average $1/2$ in the trials above).
# But numpy's implementation is using the BTPE method, which runs in $O(1)$.
# ## Conclusion
#
# - So I was write, for the inversion method the computation time is in average $O(np)$ so it is $\Omega(n)$ and cannot be $O(1)$.
# - But there has been many algorithms proposed in the literature which achieves a $O(1)$ running time, and the state-of-the-art algorithm is the BTPE method by V. Kachitvichyanukul & B. W. Schmeiser (from 1988). It is implemented in numpy, for the cases when $np > 30$ (that is, as soon as $n>60$ for $p=1/2$).
# - So the authors of [[Perturbed-History Exploration in Stochastic Multi-Armed Bandits, by Branislav Kveton, Csaba Szepesvari, Mohammad Ghavamzadeh, Craig Boutilier, 26 Feb 2019, arXiv:1902.10089]](https://arxiv.org/abs/1902.10089) were correct, it can indeed cost $O(1)$ time to generate the sum of $t$ samples from $Bern(1/2)$ (that is, a sample from $Bin(t, 1/2)$).