Section 2.3 — Inventory of discrete distributions¶

This notebook contains all the code examples from Section 2.3 Inventory of discrete distributions of the No Bullshit Guide to Statistics.

Notebook setup¶

In [1]:

# load Python modules
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# load Python modules import numpy as np import pandas as pd import seaborn as sns import matplotlib.pyplot as plt

In [2]:

# Figures setup
sns.set_theme(
    context="paper",
    style="whitegrid",
    palette="colorblind",
    rc={'figure.figsize': (6,2)},
)

%config InlineBackend.figure_format = 'retina'

# Figures setup sns.set_theme( context="paper", style="whitegrid", palette="colorblind", rc={'figure.figsize': (6,2)}, ) %config InlineBackend.figure_format = 'retina'

In [3]:

# set random seed for repeatability
np.random.seed(42)

# set random seed for repeatability np.random.seed(42)

In [4]:

%pip install -q ministats

%pip install -q ministats

[notice] A new release of pip is available: 24.3.1 -> 25.2
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.

In [5]:

from ministats import plot_pmf

from ministats import plot_pmf

Discrete uniform¶

In [6]:

from scipy.stats import randint

alpha = 1
beta = 4
rvR = randint(alpha, beta+1)

from scipy.stats import randint alpha = 1 beta = 4 rvR = randint(alpha, beta+1)

Note the +1 we added to the second argument.

In [7]:

rvR.mean()

rvR.mean()

Out[7]:

2.5

In [8]:

rvR.var()

rvR.var()

Out[8]:

1.25

The limits of the sample space of the random variable rvR can be obtained by calling its .support() method.

Probability mass function¶

In [9]:

plot_pmf(rvR, xlims=[0,8+1]);

plot_pmf(rvR, xlims=[0,8+1]);

In [10]:

rvR.support()

rvR.support()

Out[10]:

(1, 4)

Bernoulli¶

In [11]:

from scipy.stats import bernoulli

rvB = bernoulli(p=0.3)

from scipy.stats import bernoulli rvB = bernoulli(p=0.3)

In [12]:

rvB.support()

rvB.support()

Out[12]:

(0, 1)

In [13]:

rvB.mean(), rvB.var()

rvB.mean(), rvB.var()

Out[13]:

(0.3, 0.21)

In [14]:

rvB.rvs(10)

rvB.rvs(10)

Out[14]:

array([0, 1, 1, 0, 0, 0, 0, 1, 0, 1])

In [15]:

plot_pmf(rvB, xlims=[0,5]);

plot_pmf(rvB, xlims=[0,5]);

Mathematical interlude: counting combinations¶

Factorials¶

TODO explain + math formula

In [16]:

from scipy.special import factorial

factorial(4)

from scipy.special import factorial factorial(4)

Out[16]:

24.0

In [17]:

# # ALT.
# from math import factorial
# factorial(4)

# # ALT. # from math import factorial # factorial(4)

The factorial of $0$ is defined as 1: $0!=1$.

In [18]:

factorial(0)

factorial(0)

Out[18]:

1.0

The factorial function grows grows very quickly:

In [19]:

[factorial(n) for n in [5, 6, 7, 8, 9, 10, 11, 12, 13]]

[factorial(n) for n in [5, 6, 7, 8, 9, 10, 11, 12, 13]]

Out[19]:

[120.0,
 720.0,
 5040.0,
 40320.0,
 362880.0,
 3628800.0,
 39916800.0,
 479001600.0,
 6227020800.0]

Combinations¶

TODO explain + math formula

In [20]:

from scipy.special import comb

comb(5, 2)

from scipy.special import comb comb(5, 2)

Out[20]:

10.0

In [21]:

from itertools import combinations

list(combinations({1,2,3,4,5}, 2))

from itertools import combinations list(combinations({1,2,3,4,5}, 2))

Out[21]:

[(1, 2),
 (1, 3),
 (1, 4),
 (1, 5),
 (2, 3),
 (2, 4),
 (2, 5),
 (3, 4),
 (3, 5),
 (4, 5)]

Exercises¶

In [22]:

from scipy.special import comb

comb(12,6)

from scipy.special import comb comb(12,6)

Out[22]:

924.0

In [23]:

from scipy.special import comb

comb(52,5)

from scipy.special import comb comb(52,5)

Out[23]:

2598960.0

Links¶

• https://en.wikipedia.org/wiki/Factorials
• https://en.wikipedia.org/wiki/Combinations
• https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.factorial.html
• https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.comb.html
• https://docs.python.org/3/library/itertools.html#itertools.combinations

Binomial distribution¶

We'll use the name rvX because rvB was already used for the Bernoulli random variable above.

In [24]:

from scipy.stats import binom

n = 20
p = 0.14
rvX = binom(n,p)

from scipy.stats import binom n = 20 p = 0.14 rvX = binom(n,p)

In [25]:

rvX.support()

rvX.support()

Out[25]:

(0, 20)

In [26]:

rvX.mean(), rvX.var()

rvX.mean(), rvX.var()

Out[26]:

(2.8000000000000003, 2.4080000000000004)

In [27]:

plot_pmf(rvX, xlims=[0,21]);

plot_pmf(rvX, xlims=[0,21]);

Poisson¶

In [28]:

from scipy.stats import poisson
lam = 10
rvP = poisson(lam)

from scipy.stats import poisson lam = 10 rvP = poisson(lam)

In [29]:

rvP.pmf(8)

rvP.pmf(8)

Out[29]:

0.11259903214902009

In [30]:

rvP.cdf(8)

rvP.cdf(8)

Out[30]:

0.3328196787507191

In [31]:

## ALT. way to compute the value F_P(8) =
# sum([rvP.pmf(x) for x in range(0,8+1)])

## ALT. way to compute the value F_P(8) = # sum([rvP.pmf(x) for x in range(0,8+1)])

In [32]:

plot_pmf(rvP, xlims=[0,30]);

plot_pmf(rvP, xlims=[0,30]);

Geometric¶

In [33]:

from scipy.stats import geom

rvG = geom(p = 0.2)

from scipy.stats import geom rvG = geom(p = 0.2)

In [34]:

rvG.support()

rvG.support()

Out[34]:

(1, inf)

In [35]:

rvG.mean(), rvG.var()

rvG.mean(), rvG.var()

Out[35]:

(5.0, 20.0)

In [36]:

plot_pmf(rvG, xlims=[0,30]);

plot_pmf(rvG, xlims=[0,30]);

Negative binomial¶

In [37]:

from scipy.stats import nbinom

r = 10
p = 0.5
rvN = nbinom(r,p)

from scipy.stats import nbinom r = 10 p = 0.5 rvN = nbinom(r,p)

In [38]:

rvN.support()

rvN.support()

Out[38]:

(0, inf)

In [39]:

rvN.mean(), rvN.var()

rvN.mean(), rvN.var()

Out[39]:

(10.0, 20.0)

In [40]:

plot_pmf(rvN, xlims=[0,30]);

plot_pmf(rvN, xlims=[0,30]);

Computer models for discrete distributions¶

TODO table models

TODO table methods

Building computer models for random variables¶

In [41]:

# import the model family
from scipy.stats import randint

# choose parameters
alpha = 1  # start at 1
beta = 4   # stop at 4

# create the rv object
rvR = randint(alpha, beta+1)

# import the model family from scipy.stats import randint # choose parameters alpha = 1 # start at 1 beta = 4 # stop at 4 # create the rv object rvR = randint(alpha, beta+1)

Note the +1 we added to the second argument.

In [42]:

rvR.mean()

rvR.mean()

Out[42]:

2.5

In [43]:

# verify using formula
(alpha + beta) / 2

# verify using formula (alpha + beta) / 2

Out[43]:

2.5

In [44]:

rvR.var()

rvR.var()

Out[44]:

1.25

In [45]:

# verify using formula
((beta + 1 - alpha)**2 - 1) / 12

# verify using formula ((beta + 1 - alpha)**2 - 1) / 12

Out[45]:

1.25

In [46]:

rvR.std()

rvR.std()

Out[46]:

1.118033988749895

The limits of the sample space of the random variable rvR can be obtained by calling its .support() method.

In [47]:

rvR.support()

rvR.support()

Out[47]:

(1, 4)

Probability calculations¶

In [48]:

rvR.pmf(1) + rvR.pmf(2) + rvR.pmf(3)

rvR.pmf(1) + rvR.pmf(2) + rvR.pmf(3)

Out[48]:

0.75

In [49]:

[rvR.pmf(r) for r in range(1,3+1)]

[rvR.pmf(r) for r in range(1,3+1)]

Out[49]:

[0.25, 0.25, 0.25]

In [50]:

sum([rvR.pmf(r) for r in range(1,3+1)])

sum([rvR.pmf(r) for r in range(1,3+1)])

Out[50]:

0.75

In [51]:

rvR.cdf(3)

rvR.cdf(3)

Out[51]:

0.75

Plotting distributions¶

Probability mass function¶

To create a stem-plot of the probability mass function $f_R$, we can use the following three-step procedure:

1. Create a range of inputs rs for the plot.
2. Compute the value of $f_R =$ rvR for each of the inputs and store the results as list of values fRs.
3. Plot the values fRs by calling the function plt.stem(rs, fRs).

In [52]:

import matplotlib.pyplot as plt

rs = range(0, 8+1)
fRs = rvR.pmf(rs)
fRs = np.where(fRs == 0, np.nan, fRs)  # set zero fRs to np.nan
plt.stem(rs, fRs, basefmt=" ");

import matplotlib.pyplot as plt rs = range(0, 8+1) fRs = rvR.pmf(rs) fRs = np.where(fRs == 0, np.nan, fRs) # set zero fRs to np.nan plt.stem(rs, fRs, basefmt=" ");

Note, we used the np.where statement to replace zero entries in the list fRs to the special Not a Number (np.nan) values, which excludes them from the plot. This is arguably a cleaner representation of the PMF plot: the function is only defined for the sample space $\mathcal{R} = \{1,2,3,4\}$.

An alternative ,,,

In [53]:

# ALT. using helper function
from ministats import plot_pmf

plot_pmf(rvR, xlims=[0,8+1], rv_name="R");

# ALT. using helper function from ministats import plot_pmf plot_pmf(rvR, xlims=[0,8+1], rv_name="R");

Cumulative distribution function¶

In [54]:

import numpy as np
import seaborn as sns

bs = np.linspace(0, 8, 1000)
FRs = rvR.cdf(bs)
sns.lineplot(x=bs, y=FRs);

import numpy as np import seaborn as sns bs = np.linspace(0, 8, 1000) FRs = rvR.cdf(bs) sns.lineplot(x=bs, y=FRs);

In [55]:

# ALT. using helper function
from ministats import plot_cdf

plot_cdf(rvR, xlims=[0,8], rv_name="R");

# ALT. using helper function from ministats import plot_cdf plot_cdf(rvR, xlims=[0,8], rv_name="R");

Generating random observations¶

Let's generate 10 random observations from random variable rvR:

In [56]:

np.random.seed(43)

rvR.rvs(10)

np.random.seed(43) rvR.rvs(10)

Out[56]:

array([1, 1, 4, 2, 2, 3, 1, 4, 2, 4])

Explanations¶

Binomial coefficient formula explained¶

In [57]:

n = 6

[int(comb(n,k)) for k in range(0,n+1)]

n = 6 [int(comb(n,k)) for k in range(0,n+1)]

Out[57]:

[1, 6, 15, 20, 15, 6, 1]

We can draw Pascal's triangle by by repeating the above calculation for several values of n.

In [58]:

for n in range(0,10):
    row_ints = [int(comb(n,k)) for k in range(0,n+1)]
    row_str = str(row_ints)
    row_str = row_str.replace("[","").replace("]","").replace(",", "  ")
    row_str = row_str.center(60)
    print(row_str)

for n in range(0,10): row_ints = [int(comb(n,k)) for k in range(0,n+1)] row_str = str(row_ints) row_str = row_str.replace("[","").replace("]","").replace(",", " ") row_str = row_str.center(60) print(row_str)

                             1                              
                           1   1                            
                         1   2   1                          
                       1   3   3   1                        
                     1   4   6   4   1                      
                  1   5   10   10   5   1                   
                1   6   15   20   15   6   1                
             1   7   21   35   35   21   7   1              
           1   8   28   56   70   56   28   8   1           
       1   9   36   84   126   126   84   36   9   1        

Geometric sums¶

In [59]:

r = 0.3

sum([r**n for n in range(0,100)])

r = 0.3 sum([r**n for n in range(0,100)])

Out[59]:

1.4285714285714286

In [60]:

1 / (1 - r)

1 / (1 - r)

Out[60]:

1.4285714285714286

Discussion¶

Relations between distributions¶

Exercises¶

Links¶

In [ ]: