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#

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

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

%config InlineBackend.figure_format = 'retina'

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

%pip install ministats

Requirement already satisfied: ministats in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (0.3.1)
Requirement already satisfied: matplotlib>=3.7 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from ministats) (3.7.5)
Requirement already satisfied: numpy>=1.24 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from ministats) (1.24.4)
Requirement already satisfied: scipy>=1.6 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from ministats) (1.10.1)
Requirement already satisfied: pandas>=2 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from ministats) (2.0.3)
Requirement already satisfied: seaborn>=0.13.2 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from ministats) (0.13.2)
Requirement already satisfied: statsmodels>=0.14.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from ministats) (0.14.1)
Requirement already satisfied: contourpy>=1.0.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (1.1.1)
Requirement already satisfied: cycler>=0.10 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (4.51.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (24.0)
Requirement already satisfied: pillow>=6.2.0 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (10.3.0)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (2.9.0.post0)
Requirement already satisfied: importlib-resources>=3.2.0 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from matplotlib>=3.7->ministats) (6.4.0)

Requirement already satisfied: pytz>=2020.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from pandas>=2->ministats) (2024.1)
Requirement already satisfied: tzdata>=2022.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from pandas>=2->ministats) (2024.1)

Requirement already satisfied: patsy>=0.5.4 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from statsmodels>=0.14.1->ministats) (0.5.6)

Requirement already satisfied: zipp>=3.1.0 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from importlib-resources>=3.2.0->matplotlib>=3.7->ministats) (3.18.1)

Requirement already satisfied: six in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from patsy>=0.5.4->statsmodels>=0.14.1->ministats) (1.16.0)

[notice] A new release of pip is available: 23.0.1 -> 24.0
[notice] To update, run: pip install --upgrade pip

Note: you may need to restart the kernel to use updated packages.

from ministats import plot_pmf
from ministats import plot_cdf

Definitions#

Math prerequisites#

Combinatorics#

See SciPy docs:

Factorial#

from scipy.special import factorial
# ALT.
# from math import factorial

factorial(4)

24.0

 factorial(1), factorial(2), factorial(3)

(1.0, 2.0, 6.0)

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

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

factorial(15)

1307674368000.0

import numpy as np
np.log(factorial(15))/np.log(10)

12.116499611123398

Permutations#

from scipy.special import perm

perm(5,2)

20.0

perm(5,1), perm(5,2), perm(5,3), perm(5,4), perm(5,5)

(5.0, 20.0, 60.0, 120.0, 120.0)

If you want to actually see, all the possible permutations we use the permutations function from itertools

from itertools import permutations
n = 5
nitems = range(1,n+1)
k = 2
list(permutations(nitems, k))

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

Combinations#

from scipy.special import comb


comb(5,2)

10.0

from itertools import combinations

n = 5
k = 2
list(combinations(range(1,n+1), k))

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

Summations using basic Python#

N = 10
nums = range(1,N+1)
list(nums)

[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

sum(nums)

55

N*(N+1)/2  # formula

55.0

squarednums = [num**2 for num in nums]
squarednums

[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

sum(squarednums)

385

N*(N+1)*(2*N+1)/6  # formula

385.0

cubednums = [num**3 for num in nums]
cubednums

[1, 8, 27, 64, 125, 216, 343, 512, 729, 1000]

sum(cubednums)

3025

( N*(N+1)/2 )**2  # formula

3025.0

Sum of the geometric series#

a = 0.5
r = 0.5

sum([a*r**k for k in range(0,60)])

1.0

Summations using SymPy (bonus material)#

from sympy import symbols
from sympy import summation
from sympy import simplify

k, N, r = symbols("k N r")

Sum of arithmetic sequence#

a_k = k

summation(a_k, (k,0, N))
\[\displaystyle \frac{N^{2}}{2} + \frac{N}{2}\]
simplify(summation(a_k, (k,0, N)))
\[\displaystyle \frac{N \left(N + 1\right)}{2}\]

Sum of squares#

b_k = k**2

simplify(summation(b_k, (k,0, N)))
\[\displaystyle \frac{N \left(2 N^{2} + 3 N + 1\right)}{6}\]

Sum of cubes#

c_k = k**3

simplify(summation(c_k, (k,0, N)))
\[\displaystyle \frac{N^{2} \left(N^{2} + 2 N + 1\right)}{4}\]

Geometric series#

g_k = r**k

summation(g_k, (k,0, N))
\[\begin{split}\displaystyle \begin{cases} N + 1 & \text{for}\: r = 1 \\\frac{1 - r^{N + 1}}{1 - r} & \text{otherwise} \end{cases}\end{split}\]
# from sympy import limit, oo
# limit( (1-r**(N+1))/(1-r), N, oo)
# # doesn't work; need to specify assumption r < 1

Discrete distributions reference#

Discrete uniform#

# import the model family
from scipy.stats import randint

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

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

# use one of the rv object's methods

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

rvU.support()

(1, 4)

rvU.mean()

2.5

rvU.var()

1.25

rvU.std()  # = np.sqrt(rvU.var())

1.118033988749895

Probability mass function#

for x in range(1,4+1):
    print("f_U(",x,")  = ", rvU.pmf(x))

f_U( 1 )  =  0.25
f_U( 2 )  =  0.25
f_U( 3 )  =  0.25
f_U( 4 )  =  0.25

To create a stem-plot of the probability mass function \(f_U\), we can use the following three-step procedure:

  1. Create a range of inputs xs for the plot.

  2. Compute the value of \(f_U =\) rvU for each of the inputs and store the results as list of values fUs.

  3. Plot the values fUs by calling the function plt.stem(xs,fUs).

import numpy as np
import matplotlib.pyplot as plt

xs = np.arange(0,8+1)
fUs = rvU.pmf(xs)
plt.stem(xs, fUs, basefmt=" ")

<StemContainer object of 3 artists>
../_images/92ba61be193c7a280d65c4a79b05351cf11e42308e6b726c7450d89c68c2eac6.png 
# ALT
plot_pmf(rvU, xlims=[0,8+1])

<Axes: xlabel='x', ylabel='$f_{X}$'>
../_images/86e80c9cec5d40d6e52dc02ffac854d52a4ed793935bf553a6f848485d4d6945.png

Cumulative distribution function#

for b in range(1,4+1):
    print("F_U(", b, ")  = ", rvU.cdf(b))

F_U( 1 )  =  0.25
F_U( 2 )  =  0.5
F_U( 3 )  =  0.75
F_U( 4 )  =  1.0

import numpy as np
import seaborn as sns

xs = np.linspace(0,8,1000)
FUs = rvU.cdf(xs)
sns.lineplot(x=xs, y=FUs)

<Axes: >
../_images/414377e02499ded99c63382044ccc2107b89948d4593b2634348c038e3c1c875.png 
# ALT
plot_cdf(rvU, xlims=[0,8], rv_name="U")

<Axes: xlabel='u', ylabel='$F_{U}$'>
../_images/e7476c895474c00a26b5f47a290af6b0dbd2bd50ee9476ada54ff32b386229b5.png

Let’s generate 10 random observations from random variable rvU:

rvU.rvs(10)

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

Bernoulli#

from scipy.stats import bernoulli

rvB = bernoulli(p=0.3)

rvB.support()

(0, 1)

rvB.mean(), rvB.var()

(0.3, 0.21)

rvB.rvs(10)

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

plot_pmf(rvB, xlims=[0,5]);
../_images/e12834fae913001510db8464ff1e9bb75bfb33ef5f9cbfa92af41c71bc666ab2.png

Poisson#

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

rvP.pmf(8)

0.11259903214902009

rvP.cdf(8)

0.3328196787507191

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

plot_pmf(rvP, xlims=[0,30]);
../_images/539964ad94640fae2985970ab788a898243d74e881ad8aa0f5c956f33425ea36.png

Binomial#

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

from scipy.stats import binom

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

rvX.support()

(0, 20)

rvX.mean(), rvX.var()

(2.8000000000000003, 2.4080000000000004)

plot_pmf(rvX, xlims=[0,30]);
../_images/eeac252813578c6428cef501b9156df4e32c99a4b42cc25e40b4dae5cd42b8da.png

Geometric#

from scipy.stats import geom

rvG = geom(p = 0.2)

rvG.support()

(1, inf)

rvG.mean(), rvG.var()

(5.0, 20.0)

plot_pmf(rvG, xlims=[0,40]);
../_images/59cb300b61db6111fa927118f70790ed646fdcd6f2aae1214063b705e03816ba.png

Negative binomial#

from scipy.stats import nbinom

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

rvN.support()

(0, inf)

rvN.mean(), rvN.var()

(10.0, 20.0)

plot_pmf(rvN, xlims=[0,40]);
../_images/41269d859ca2432e2d05bd2bc5ff8f454ac1c417a38fcc6d51f4037cbd2a76e7.png

Hypergeometric#

from scipy.stats import hypergeom

a = 30   # number of success balls
b = 40   # number of failure balls
n = 20   # how many we're drawing

rvH = hypergeom(a+b, a, n)

rvH.support()

(0, 20)

rvH.mean(), rvH.var()

(8.571428571428571, 3.54924578527063)

meanH, stdH = rvH.stats()
print("mean =", meanH, "  std =", stdH)

mean = 8.571428571428571   std = 3.54924578527063

plot_pmf(rvH, xlims=[0,30]);
../_images/f8183ec7137c5070902ae526315910b50112fad7dbdebf49c53dfa0955bed3e6.png

Tomatoes salad probabilities#

a = 3   # number of good tomatoes
b = 4   # number of rotten tomatoes
n = 2   # how many we're drawing

rvHe = hypergeom(a+b, a, n)


plot_pmf(rvHe, xlims=[0,3])

rvHe.pmf(0), rvHe.pmf(1), rvHe.pmf(2)

(0.28571428571428575, 0.5714285714285715, 0.14285714285714288)
../_images/a958763b0f96047f115f579c126e2978514c93da2a830f2571b9546d146c15d7.png

Number of dogs seen by Amy#

a = 7        # number dogs
b = 20 - 7   # number of other animals
n = 12       # how many "patients" Amy will see today

rvD = hypergeom(a+b, a, n)

# Pr of exactly five dogs
rvD.pmf(5)

0.2860681114551084

plot_pmf(rvD, xlims=[0,10]);
../_images/b3cd85f41a4c1fea96f76defea4d70ea25c5b92318bd65e82a532be801c5bf83.png

Multinomial#

See docs.

from scipy.stats import multinomial

n = 10
ps = [0.1, 0.5, 0.8]

rvM = multinomial(n,ps)

rvM.rvs()

array([[0, 6, 4]])

# TODO: 3D scatter plot of points in space

Modelling real-world data using probability#

TODO: add simple inference and plots

Discussion#