Univariate Discrete Distributions¶
Statistical Distributions — the discrete family, from scratch¶
Folder 1 of the Statistical Distributions catalog. For each distribution we implement the probability mass function (PMF), cumulative distribution (CDF), quantile, random-number generator, and moments from scratch — using only elementary math and the (log-)gamma function — and validate against scipy.stats as an independent reference. The R companion validates against R's d/p/q/r functions.
What defines a discrete distribution¶
A discrete distribution places probability mass on a countable set (its support). Three objects describe it:
- the PMF $p(k)=\Pr(X=k)$, with $\sum_k p(k)=1$;
- the CDF $F(k)=\Pr(X\le k)=\sum_{j\le k}p(j)$, a right-continuous step function;
- the moments — mean $\mu=\mathbb E[X]$, variance $\sigma^2$, skewness (asymmetry) and excess kurtosis (tail-heaviness relative to normal), obtainable from the probability generating function $G(z)=\mathbb E[z^X]$.
Every distribution below is checked three ways: from-scratch vs scipy (PMF & CDF agree to machine precision), analytical vs empirical moments (closed form vs the mean/var/skew/kurtosis of a large sample), and RNG vs PMF (a histogram of simulated draws against the analytical mass). The random-number-generation algorithm is itself part of the story — inverse-transform, Bernoulli sums, Knuth's Poisson method, compound mixtures.
Families covered¶
- Foundations — Bernoulli, Discrete Uniform, Categorical
- The counting family — Binomial, Poisson (and the Binomial→Poisson limit)
- Waiting times — Geometric, Negative Binomial
- Sampling without replacement — Hypergeometric
- Over-dispersion & flexible dispersion — Beta-Binomial, Zero-Inflated Poisson, Conway–Maxwell–Poisson
- Differences & heavy tails — Skellam, Logarithmic, Yule–Simon
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import discrete as D
from scipy import stats as S
plt.rcParams.update({"figure.figsize": (9,4), "axes.grid": True, "grid.alpha": .22,
"axes.spines.top": False, "axes.spines.right": False, "font.size": 10.5})
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#718096","#6b46c1"
rng = np.random.default_rng(0)
def emp_moments(x):
x = np.asarray(x, float); m = x.mean(); c2 = x.var()
return m, c2, ((x-m)**3).mean()/c2**1.5, ((x-m)**4).mean()/c2**2 - 3
def panel(name, ks, pmf_fs, cdf_fs, samples, m_fs, sp_pmf=None, sp_cdf=None, m_sp=None):
'''Standard validation figure + moments table for one discrete distribution.'''
fig, ax = plt.subplots(1, 3, figsize=(14.5, 3.5))
ax[0].vlines(ks, 0, pmf_fs, color=BLUE, lw=2, alpha=.8)
if sp_pmf is not None: ax[0].plot(ks, sp_pmf, "o", color=RED, ms=4.5, mfc="none", label="scipy")
ax[0].plot(ks, pmf_fs, ".", color=BLUE, ms=3, label="from-scratch"); ax[0].set_title(name+" · PMF")
ax[0].legend(fontsize=8); ax[0].set_xlabel("k")
ax[1].step(np.append(ks, ks[-1]+1), np.append(cdf_fs, cdf_fs[-1]), where="post", color=BLUE, lw=2)
if sp_cdf is not None: ax[1].plot(ks, sp_cdf, ".", color=RED, ms=6)
ax[1].set_title("CDF (step = from-scratch, dots = scipy)"); ax[1].set_ylim(-.03,1.03); ax[1].set_xlabel("k")
edges = np.arange(ks.min()-.5, ks.max()+1.5)
ax[2].hist(samples, bins=edges, density=True, color="#cfe3f6", edgecolor="white", lw=.4, label="RNG samples")
ax[2].plot(ks, pmf_fs, "o-", color=BLUE, ms=3, lw=1, label="analytical PMF")
ax[2].set_title("RNG vs analytical (N=%s)" % f"{len(samples):,}"); ax[2].legend(fontsize=8); ax[2].set_xlabel("k")
plt.tight_layout(); plt.show()
me = emp_moments(samples)
cols = {"analytical (formula)": [m_fs["mean"], m_fs["var"], m_fs["skew"], m_fs["exkurt"]]}
if m_sp is not None: cols["scipy"] = list(m_sp)
cols["empirical (RNG)"] = list(me)
print(pd.DataFrame(cols, index=["mean","variance","skewness","excess kurtosis"]).round(3).to_string())
def sp_stats(dist): return [float(x) for x in dist.stats("mvsk")]
print("helpers ready")
helpers ready
1. Foundations¶
Bernoulli(p) — a single yes/no trial, with $p(1)=p$ and $p(0)=1-p$. The atom from which the Binomial is built; mean $p$, variance $p(1-p)$.
Discrete Uniform{a,…,b} — every integer in a range equally likely; the maximum-entropy discrete law on a bounded support; symmetric (skew $0$), platykurtic.
Categorical(π₁,…,π_K) — one draw over $K$ labelled outcomes; the atom of the Multinomial and the likelihood behind classification.
# Bernoulli(0.35)
ks = np.array([0,1]); p = 0.35
panel("Bernoulli(0.35)", ks, D.bernoulli_pmf(ks,p), D.bernoulli_cdf(ks,p),
D.bernoulli_rng(p,100000,rng), D.bernoulli_moments(p),
S.bernoulli.pmf(ks,p), S.bernoulli.cdf(ks,p), sp_stats(S.bernoulli(p)))
analytical (formula) scipy empirical (RNG) mean 0.350 0.350 0.348 variance 0.227 0.227 0.227 skewness 0.629 0.629 0.639 excess kurtosis -1.604 -1.604 -1.591
# Discrete Uniform {2,...,8}
a,b = 2,8; ks = np.arange(a,b+1)
panel("Discrete Uniform{2..8}", ks, D.duniform_pmf(ks,a,b), D.duniform_cdf(ks,a,b),
D.duniform_rng(a,b,100000,rng), D.duniform_moments(a,b),
S.randint.pmf(ks,a,b+1), S.randint.cdf(ks,a,b+1), sp_stats(S.randint(a,b+1)))
analytical (formula) scipy empirical (RNG) mean 5.00 5.00 4.994 variance 4.00 4.00 3.993 skewness 0.00 0.00 0.005 excess kurtosis -1.25 -1.25 -1.246
# Categorical over 5 classes
probs = np.array([.30,.10,.25,.20,.15]); ks = np.arange(len(probs))
m = D.moments_from_pmf(ks, D.categorical_pmf(ks,probs))
panel("Categorical (K=5)", ks, D.categorical_pmf(ks,probs), D.categorical_cdf(ks,probs),
D.categorical_rng(probs,100000,rng), m)
print("\n(No scipy Categorical; PMF is the defining probability vector, moments by summation.)")
analytical (formula) empirical (RNG) mean 1.800 1.801 variance 2.060 2.061 skewness 0.049 0.046 excess kurtosis -1.322 -1.323 (No scipy Categorical; PMF is the defining probability vector, moments by summation.)
2. The counting family: Binomial & Poisson¶
Binomial(n, p) counts successes in $n$ independent Bernoulli trials: $p(k)=\binom{n}{k}p^k(1-p)^{n-k}$. Mean $np$, variance $np(1-p)$ (always under-dispersed relative to Poisson).
Poisson(λ) counts events in a fixed interval when they occur independently at constant rate: $p(k)=e^{-\lambda}\lambda^k/k!$. Its defining feature is equidispersion, mean $=$ variance $=\lambda$; skewness $1/\sqrt\lambda$ and excess kurtosis $1/\lambda$ both vanish as $\lambda\to\infty$ (it approaches Normal).
The two are linked by a limit: $\mathrm{Binomial}(n,\lambda/n)\to\mathrm{Poisson}(\lambda)$ as $n\to\infty$ — rare events among many trials. Our RNGs mirror the definitions: Binomial as a sum of Bernoullis, Poisson by Knuth's algorithm (multiply uniforms until the product drops below $e^{-\lambda}$).
n,p = 20,0.35; ks = np.arange(0,n+1)
panel("Binomial(20, 0.35)", ks, D.binomial_pmf(ks,n,p), D.binomial_cdf(ks,n,p),
D.binomial_rng(n,p,150000,rng), D.binomial_moments(n,p),
S.binom.pmf(ks,n,p), S.binom.cdf(ks,n,p), sp_stats(S.binom(n,p)))
analytical (formula) scipy empirical (RNG) mean 7.000 7.000 6.999 variance 4.550 4.550 4.532 skewness 0.141 0.141 0.139 excess kurtosis -0.080 -0.080 -0.073
lam = 4.0; ks = np.arange(0,16)
panel("Poisson(4)", ks, D.poisson_pmf(ks,lam), D.poisson_cdf(ks,lam),
D.poisson_rng(lam,100000,rng), D.poisson_moments(lam),
S.poisson.pmf(ks,lam), S.poisson.cdf(ks,lam), sp_stats(S.poisson(lam)))
analytical (formula) scipy empirical (RNG) mean 4.00 4.00 3.999 variance 4.00 4.00 3.998 skewness 0.50 0.50 0.502 excess kurtosis 0.25 0.25 0.279
# Feature plots: Binomial -> Poisson limit, and Poisson's central-limit drift
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
lam = 3.0; kk = np.arange(0,14)
for n_,c in zip([6,12,30,120], [ORANGE,GREEN,PURP,BLUE]):
ax[0].plot(kk, D.binomial_pmf(kk, n_, lam/n_), "o-", ms=3, lw=1, color=c, alpha=.8, label="Binomial(n=%d)"%n_)
ax[0].plot(kk, D.poisson_pmf(kk, lam), "k--", lw=2, label="Poisson(3) limit")
ax[0].set_title("Binomial(n, 3/n) → Poisson(3) as n grows"); ax[0].set_xlabel("k"); ax[0].legend(fontsize=8)
lams = np.linspace(0.5, 30, 60)
ax[1].plot(lams, 1/np.sqrt(lams), color=BLUE, lw=2, label="skewness $1/\\sqrt{\\lambda}$")
ax[1].plot(lams, 1/lams, color=RED, lw=2, label="excess kurtosis $1/\\lambda$")
ax[1].axhline(0, color=GREY, lw=.8); ax[1].set_xlabel("Poisson rate $\\lambda$"); ax[1].set_ylabel("moment")
ax[1].set_title("Poisson approaches Normal as $\\lambda$ grows"); ax[1].legend()
plt.tight_layout(); plt.show()
3. Waiting times: Geometric & Negative Binomial¶
Geometric(p) — the number of trials to the first success, $p(k)=(1-p)^{k-1}p$, $k=1,2,\dots$ The unique memoryless discrete law; heavy geometric (exponential) tail. Sampled in closed form by inverse transform, $k=\lceil\ln U/\ln(1-p)\rceil$.
Negative Binomial(r, p) — the number of failures before the $r$-th success, $p(k)=\binom{k+r-1}{k}p^r(1-p)^k$. Crucially it is also the Poisson–Gamma mixture: draw a rate $\lambda\sim\text{Gamma}(r,\cdot)$, then $k\sim\text{Poisson}(\lambda)$. That mixing makes it over-dispersed (variance $>$ mean) — the standard fix for count data too variable for Poisson. Our RNG uses exactly this mixture.
p = 0.3; ks = np.arange(1,25)
panel("Geometric(0.3)", ks, D.geometric_pmf(ks,p), D.geometric_cdf(ks,p),
D.geometric_rng(p,100000,rng), D.geometric_moments(p),
S.geom.pmf(ks,p), S.geom.cdf(ks,p), sp_stats(S.geom(p)))
analytical (formula) scipy empirical (RNG) mean 3.333 3.333 3.326 variance 7.778 7.778 7.689 skewness 2.032 2.032 1.996 excess kurtosis 6.129 6.129 5.766
r,p = 5,0.4; ks = np.arange(0,35)
panel("NegBinomial(r=5, p=0.4)", ks, D.nbinom_pmf(ks,r,p), D.nbinom_cdf(ks,r,p),
D.nbinom_rng(r,p,120000,rng), D.nbinom_moments(r,p),
S.nbinom.pmf(ks,r,p), S.nbinom.cdf(ks,r,p), sp_stats(S.nbinom(r,p)))
analytical (formula) scipy empirical (RNG) mean 7.500 7.500 7.509 variance 18.750 18.750 18.882 skewness 0.924 0.924 0.916 excess kurtosis 1.253 1.253 1.208
# Feature: over-dispersion -- Poisson vs Negative Binomial at the SAME mean
mean = 7.5; kk = np.arange(0,30) # NegBinom(r=5, p=5/(5+mean)) has mean r(1-p)/p = 7.5
plt.figure(figsize=(9,4))
plt.vlines(kk-0.12, 0, D.poisson_pmf(kk, mean), color=BLUE, lw=2.5, label="Poisson(7.5) var=7.5")
plt.vlines(kk+0.12, 0, D.nbinom_pmf(kk, 5, 5/(5+mean)), color=RED, lw=2.5, label="NegBinom mean=7.5, var=%.1f"%D.nbinom_moments(5,5/(5+mean))["var"])
plt.title("Same mean, different spread: the Negative Binomial is over-dispersed")
plt.xlabel("k"); plt.legend(); plt.show()
print("Equal means (7.5); the Poisson-Gamma mixing inflates the NegBinom variance and fattens its right tail.")
Equal means (7.5); the Poisson-Gamma mixing inflates the NegBinom variance and fattens its right tail.
4. Sampling without replacement: Hypergeometric¶
Hypergeometric(M, K, n) — draw $n$ items without replacement from a population of $M$ containing $K$ successes; count the successes: $p(k)=\binom{K}{k}\binom{M-K}{n-k}/\binom{M}{n}$. It is the without-replacement analogue of the Binomial; the finite population makes draws dependent, so its variance carries the finite-population correction $\frac{M-n}{M-1}<1$ — sampling without replacement is less variable. As $M\to\infty$ it converges to the Binomial. Our RNG shuffles an explicit urn.
M,K,n = 50,15,10; ks = np.arange(0,min(n,K)+1)
panel("Hypergeometric(M=50,K=15,n=10)", ks, D.hypergeom_pmf(ks,M,K,n), D.hypergeom_cdf(ks,M,K,n),
D.hypergeom_rng(M,K,n,25000,rng), D.hypergeom_moments(M,K,n),
S.hypergeom.pmf(ks,M,K,n), S.hypergeom.cdf(ks,M,K,n), sp_stats(S.hypergeom(M,K,n)))
# contrast with Binomial (with replacement) at the same success fraction
kk = np.arange(0,11)
plt.figure(figsize=(9,3.8))
plt.vlines(kk-0.12, 0, D.hypergeom_pmf(kk,M,K,n), color=BLUE, lw=2.5, label="Hypergeometric (no replacement)")
plt.vlines(kk+0.12, 0, D.binomial_pmf(kk,n,K/M), color=RED, lw=2.5, label="Binomial(10, 0.3) (with replacement)")
plt.title("Without vs with replacement — hypergeometric is tighter (finite-population correction)")
plt.xlabel("k"); plt.legend(); plt.show()
analytical (formula) scipy empirical (RNG) mean 3.000 3.000 3.002 variance 1.714 1.714 1.731 skewness 0.191 0.191 0.198 excess kurtosis -0.127 -0.127 -0.116
5. Over-dispersion & flexible dispersion¶
Beta-Binomial(n, a, b) — a Binomial whose success probability is itself random, $p\sim\text{Beta}(a,b)$: $p(k)=\binom{n}{k}\frac{B(k+a,\,n-k+b)}{B(a,b)}$. The random $p$ adds variance — the over-dispersed count model for bounded counts.
Zero-Inflated Poisson(π, λ) — a mixture that inflates the zero: with probability $\pi$ a structural zero, otherwise a Poisson draw. Models count data with more zeros than Poisson allows.
Conway–Maxwell–Poisson(λ, ν) — a two-parameter generalisation with a dispersion knob $\nu$: $\nu=1$ recovers Poisson, $\nu<1$ gives over-dispersion, $\nu>1$ under-dispersion. The normalising constant $Z(\lambda,\nu)=\sum_k \lambda^k/(k!)^\nu$ has no closed form, so we compute it by series — a good example of a distribution defined only up to a computable constant.
n,a,b = 20,2,5; ks = np.arange(0,n+1)
panel("Beta-Binomial(20, a=2, b=5)", ks, D.betabinom_pmf(ks,n,a,b), D.betabinom_cdf(ks,n,a,b),
D.betabinom_rng(n,a,b,100000,rng), D.betabinom_moments(n,a,b),
S.betabinom.pmf(ks,n,a,b), S.betabinom.cdf(ks,n,a,b), sp_stats(S.betabinom(n,a,b)))
analytical (formula) scipy empirical (RNG) mean 5.714 5.714 5.715 variance 13.776 13.776 13.824 skewness 0.603 0.603 0.598 excess kurtosis -0.135 -0.135 -0.146
pi,lam = 0.3,4.0; ks = np.arange(0,16)
ref_pmf = np.where(ks==0, pi + (1-pi)*np.exp(-lam), (1-pi)*S.poisson.pmf(ks,lam)) # explicit mixture reference
ref_cdf = pi + (1-pi)*S.poisson.cdf(ks,lam)
panel("Zero-Inflated Poisson(π=0.3, λ=4)", ks, D.zip_pmf(ks,pi,lam), D.zip_cdf(ks,pi,lam),
D.zip_rng(pi,lam,100000,rng), D.zip_moments(pi,lam), ref_pmf, ref_cdf)
print("\n(Reference = explicit π·δ₀ + (1-π)·Poisson mixture; note the mass spike at k=0.)")
analytical (formula) empirical (RNG) mean 2.800 2.801 variance 6.160 6.142 skewness 0.491 0.486 excess kurtosis -0.525 -0.529 (Reference = explicit π·δ₀ + (1-π)·Poisson mixture; note the mass spike at k=0.)
# COM-Poisson: dispersion knob nu
ks = np.arange(0,20)
panel("COM-Poisson(λ=6, ν=1.5, under-disp.)", ks, D.compoisson_pmf(ks,6.,1.5), D.compoisson_cdf(ks,6.,1.5),
D.compoisson_rng(6.,1.5,80000,rng), D.compoisson_moments(6.,1.5))
kk = np.arange(0,22)
plt.figure(figsize=(9,4))
for nu,c,lab in [(0.5,RED,"ν=0.5 (over-disp.)"),(1.0,GREY,"ν=1 (Poisson)"),(1.8,BLUE,"ν=1.8 (under-disp.)")]:
plt.plot(kk, D.compoisson_pmf(kk,6.,nu), "o-", ms=3, lw=1, color=c, label=lab)
plt.title("Conway–Maxwell–Poisson: ν tunes dispersion around Poisson (λ=6)"); plt.xlabel("k"); plt.legend(); plt.show()
analytical (formula) empirical (RNG) mean 3.126 3.132 variance 2.210 2.230 skewness 0.444 0.453 excess kurtosis 0.202 0.246
6. Differences & heavy tails¶
Skellam(μ₁, μ₂) — the difference of two independent Poissons (e.g. a sports score margin); supported on all integers, mean $\mu_1-\mu_2$, variance $\mu_1+\mu_2$. Its PMF involves a modified Bessel function.
Logarithmic / log-series (p) — $p(k)\propto p^k/k$, $k\ge1$; the classic species-abundance law and the zero-truncated limit of the Negative Binomial.
Yule–Simon(ρ) — a power-law distribution, $p(k)=\rho\,B(k,\rho+1)\sim k^{-(\rho+1)}$. Its heavy polynomial tail (vs the geometric's exponential tail) models word frequencies, city sizes, citations — the "rich get richer" processes. It has finite mean only for $\rho>1$, finite variance only for $\rho>2$.
mu1,mu2 = 3.0,2.0; ks = np.arange(-12,15)
panel("Skellam(μ₁=3, μ₂=2)", ks, D.skellam_pmf(ks,mu1,mu2), D.skellam_cdf(ks,mu1,mu2),
D.skellam_rng(mu1,mu2,100000,rng), D.skellam_moments(mu1,mu2),
S.skellam.pmf(ks,mu1,mu2), S.skellam.cdf(ks,mu1,mu2), sp_stats(S.skellam(mu1,mu2)))
analytical (formula) scipy empirical (RNG) mean 1.000 1.000 1.001 variance 5.000 5.000 5.013 skewness 0.089 0.089 0.091 excess kurtosis 0.200 0.200 0.201
p = 0.6; ks = np.arange(1,20)
panel("Logarithmic(0.6)", ks, D.logarithmic_pmf(ks,p), D.logarithmic_cdf(ks,p),
D.logarithmic_rng(p,100000,rng), D.logarithmic_moments(p),
S.logser.pmf(ks,p), S.logser.cdf(ks,p), sp_stats(S.logser(p)))
analytical (formula) scipy empirical (RNG) mean 1.637 1.637 1.637 variance 1.413 1.413 1.383 skewness 3.005 3.005 2.907 excess kurtosis 13.656 13.656 12.472
rho = 2.5; ks = np.arange(1,25)
panel("Yule–Simon(ρ=2.5)", ks, D.yule_pmf(ks,rho), D.yule_cdf(ks,rho),
D.yule_rng(rho,100000,rng), D.yule_moments(rho),
S.yulesimon.pmf(ks,rho), S.yulesimon.cdf(ks,rho), sp_stats(S.yulesimon(rho)))
# Feature: power-law vs exponential tail on log-log axes
kk = np.arange(1,400)
plt.figure(figsize=(9,4))
plt.loglog(kk, D.yule_pmf(kk,2.5), color=BLUE, lw=2, label="Yule–Simon (power-law tail)")
plt.loglog(kk, D.geometric_pmf(kk,0.15), color=RED, lw=2, label="Geometric (exponential tail)")
plt.xlabel("k (log)"); plt.ylabel("P(X=k) (log)"); plt.title("Heavy vs light tail: Yule–Simon is a straight line on log-log")
plt.legend(); plt.show()
analytical (formula) scipy empirical (RNG) mean 1.667 1.667 1.666 variance 5.556 5.556 4.859 skewness inf inf 22.860 excess kurtosis inf inf 1336.455
7. Timing: the price of transparency¶
For a fixed number of draws, which approach is faster? We time both, and the answer splits by algorithm. Where our method is itself vectorised, lean numpy often beats scipy.stats's general-purpose .rvs — which pays overhead for its generic machinery. Where our method uses an explicit Python loop, it is 1–2 orders of magnitude slower than the library's compiled sampler. Correctness is identical throughout; only speed differs.
import time
def best_time(fn, reps=3):
t = np.inf
for _ in range(reps):
s = time.perf_counter(); fn(); t = min(t, time.perf_counter()-s)
return t
N = 100000
benches = [
("Bernoulli", lambda: D.bernoulli_rng(.35,N,rng), lambda: S.bernoulli.rvs(.35,size=N,random_state=rng), "vectorised"),
("Geometric", lambda: D.geometric_rng(.3,N,rng), lambda: S.geom.rvs(.3,size=N,random_state=rng), "vectorised (inverse-CDF)"),
("Binomial", lambda: D.binomial_rng(20,.35,N,rng), lambda: S.binom.rvs(20,.35,size=N,random_state=rng), "vectorised (Bernoulli sum)"),
("NegBinom", lambda: D.nbinom_rng(5,.4,N,rng), lambda: S.nbinom.rvs(5,.4,size=N,random_state=rng), "vectorised (Poisson–Gamma)"),
("Beta-Binom", lambda: D.betabinom_rng(20,2,5,N,rng), lambda: S.betabinom.rvs(20,2,5,size=N,random_state=rng), "vectorised (compound)"),
("Skellam", lambda: D.skellam_rng(3,2,N,rng), lambda: S.skellam.rvs(3,2,size=N,random_state=rng), "vectorised (Poisson diff.)"),
("Yule–Simon", lambda: D.yule_rng(2.5,N,rng), lambda: S.yulesimon.rvs(2.5,size=N,random_state=rng), "vectorised (exp–geom mix)"),
("Poisson", lambda: D.poisson_rng(4.,N,rng), lambda: S.poisson.rvs(4.,size=N,random_state=rng), "PYTHON LOOP (Knuth)"),
("Hypergeom", lambda: D.hypergeom_rng(50,15,10,N//5,rng),lambda: S.hypergeom.rvs(50,15,10,size=N//5,random_state=rng),"PYTHON LOOP (urn)"),
]
rows = []
for name, fs, sp, kind in benches:
nn = N//5 if "Hypergeom" in name else N
tf, ts = best_time(fs), best_time(sp)
rows.append((name, nn/tf/1e6, nn/ts/1e6, ts/tf, kind))
df = pd.DataFrame(rows, columns=["dist","from-scratch (M/s)","scipy (M/s)","ratio","method"]).set_index("dist")
print(df.round(2).to_string())
order = df["from-scratch (M/s)"].index
y = np.arange(len(order)); h = .38
loop = df["method"].str.contains("LOOP")
fig, ax = plt.subplots(figsize=(11,4.6))
ax.barh(y+h/2, df["from-scratch (M/s)"], h, color=np.where(loop, RED, BLUE), label="from-scratch")
ax.barh(y-h/2, df["scipy (M/s)"], h, color=GREY, alpha=.8, label="scipy / numpy")
ax.set_yticks(y); ax.set_yticklabels(order); ax.set_xscale("log")
ax.set_xlabel("throughput — million samples / second (log)"); ax.invert_yaxis()
ax.set_title("RNG throughput: vectorised from-scratch (blue) keeps up; Python-loop methods (red) do not")
ax.legend(); plt.tight_layout(); plt.show()
# summary text computed from the table above -- never hardcode a timing claim
wins = df[(df["ratio"] > 1.15) & ~df["method"].str.contains("LOOP")]["ratio"]
loops = df[df["method"].str.contains("LOOP")]["ratio"]
print("\nFor a fixed N, the vectorised from-scratch generators that BEAT scipy's generic .rvs: "
+ ", ".join(f"{d} ~{r:.1f}x" for d, r in wins.items()) + ".")
print("Lean numpy avoids scipy's wrapper overhead; the remaining vectorised methods sit at parity.")
print(f"Only the explicit Python loops (Knuth's Poisson, urn Hypergeometric) lose badly "
f"({1/loops.max():.0f}-{1/loops.min():.0f}x slower):")
print("that is the true cost of one-draw-at-a-time transparency, and why libraries compile their samplers.")
from-scratch (M/s) scipy (M/s) ratio method dist Bernoulli 558.04 106.83 5.22 vectorised Geometric 195.27 107.17 1.82 vectorised (inverse-CDF) Binomial 13.30 32.56 0.41 vectorised (Bernoulli sum) NegBinom 19.98 19.17 1.04 vectorised (Poisson–Gamma) Beta-Binom 13.18 13.03 1.01 vectorised (compound) Skellam 22.01 20.20 1.09 vectorised (Poisson diff.) Yule–Simon 71.22 60.26 1.18 vectorised (exp–geom mix) Poisson 0.80 36.09 0.02 PYTHON LOOP (Knuth) Hypergeom 0.38 21.26 0.02 PYTHON LOOP (urn)
For a fixed N, the vectorised from-scratch generators that BEAT scipy's generic .rvs: Bernoulli ~5.2x, Geometric ~1.8x, Yule–Simon ~1.2x. Lean numpy avoids scipy's wrapper overhead; the remaining vectorised methods sit at parity. Only the explicit Python loops (Knuth's Poisson, urn Hypergeometric) lose badly (45-56x slower): that is the true cost of one-draw-at-a-time transparency, and why libraries compile their samplers.
8. Summary¶
Fourteen discrete distributions, each with a from-scratch PMF / CDF / quantile / RNG / moments that matches scipy.stats to machine precision, closed-form moments that match a large simulated sample, and a random-number generator whose algorithm reflects the distribution's structure (Bernoulli sums, Knuth's Poisson, Poisson–Gamma mixing, urn shuffling, inverse transform).
The relationships are the point. Bernoulli builds the Binomial; the Binomial's rare-event limit is the Poisson; mixing a Poisson rate over a Gamma gives the over-dispersed Negative Binomial (and truncating it, the Logarithmic); a Beta-random probability gives the Beta-Binomial; drawing without replacement gives the Hypergeometric; the difference of Poissons gives the Skellam; and a preferential-attachment process gives the power-law Yule–Simon. Conway–Maxwell–Poisson sits above the Poisson with an explicit dispersion dial.
Next in the catalog: the continuous families — beginning with the Gaussian & exponential family, then heavy tails, extreme value / survival, and the Bayesian prior & augmentation toolkit.