Statistical Distributions
A catalog built on one rule: implement everything from scratch. For each distribution the density or mass function, the CDF, the quantile, the random-number generator and the moments are written from elementary math and the log-gamma function — then validated against an independent reference library, and re-implemented a second time in another language against a second reference. Agreement across two languages and two libraries is what separates a correct formula from a merely self-consistent one. The sampling algorithms are treated as content rather than plumbing — inverse transform, compound mixtures, Knuth's Poisson, urn shuffling — and each notebook closes by timing them against the compiled library samplers, which prices the cost of doing it transparently. The catalog begins with the discrete families; the continuous ones — the Gaussian & exponential family, heavy tails & skew, extreme value & survival, and the Bayesian prior & augmentation toolkit — follow.
Univariate Discrete Distributions
Fourteen discrete laws — Bernoulli, Discrete Uniform, Categorical, Binomial,
Poisson, Geometric, Negative Binomial, Hypergeometric, Beta-Binomial, Zero-Inflated
Poisson, Conway–Maxwell–Poisson, Skellam, Logarithmic and Yule–Simon
— each with its PMF, CDF, quantile, RNG and moments written from scratch, validated against
scipy.stats in Python and independently against base R's d/p/q/r
plus extraDistr in R. The organising theme is how the families
generate one another: Bernoulli builds the Binomial, whose rare-event limit is the
Poisson; mixing a Poisson rate over a Gamma gives the over-dispersed Negative Binomial, and
truncating that gives the Logarithmic; a Beta-random success probability gives the
Beta-Binomial; sampling without replacement gives the Hypergeometric and its
finite-population correction; differencing two Poissons gives the Skellam; and preferential
attachment gives the power-law Yule–Simon, whose skewness and kurtosis are genuinely
infinite at the illustrated ρ. Each panel carries a triple check — from-scratch versus
library, analytical versus empirical moments, and simulated draws versus the analytical mass —
and both notebooks close by benchmarking RNG throughput, where vectorised from-scratch code
runs at parity with the compiled samplers while the explicit per-draw loops (Knuth's Poisson,
the urn Hypergeometric) run one to two orders of magnitude slower.
The Gaussian & Exponential Family
The continuous counterpart to the discrete catalog — ten laws spanning the
Uniform, Normal, Exponential, Gamma, Inverse-Gamma, Beta, Chi-square,
Student-t, F and Log-normal, each with its PDF, CDF, RNG and
moments from scratch, validated against scipy.stats and independently in base R
— where every CDF is recovered by integrating the notebook's own PDF rather than
calling pgamma or pbeta, so the check is not circular. The spine is
a construction chain in which no generator calls a library: two uniforms become a Normal
(Box–Muller), a uniform becomes an Exponential, Normal plus uniform build the Gamma
(Marsaglia–Tsang rejection), and the Gamma then spawns the Chi-square, the Inverse-Gamma and
— as a ratio — the Beta, from which Student-t, F and the Log-normal follow. A dedicated
section writes seven of the ten in exponential-family canonical form
(natural parameter, sufficient statistic, base measure, log-partition), which is exactly why
those carry conjugate priors — while the Uniform, t and F sit outside it. Throughout, moments
that do not exist are reported as undefined rather than approximated: the Inverse-Gamma at
a=3 has a finite mean and variance but no skewness or kurtosis.
Heavy Tails & Skewness
The laws built for what the Gaussian misses — Laplace, Logistic, Cauchy, Pareto, Lévy, Student-t, the Generalized Error Distribution, the Skew-Normal and Hansen's skew-t, plus the α-stable family that unifies them. The organising fact is that heavy enough tails destroy moments, and the notebooks treat that as content rather than an inconvenience: the Cauchy gets no moment table at all because it has none, the Lévy's mean is infinite, and a Pareto reports its skewness as undefined because it exists only below the tail index. The centrepiece is the α-stable law, which has no closed-form density in general — it is sampled by the Chambers–Mallows–Stuck transform and validated through its characteristic function rather than a PDF, with one sampler shown reproducing the Normal, Cauchy and Lévy as its α = 2, 1 and ½ cases. Hansen's skew-t closes it as a finance capstone, standardized to mean 0 and variance 1 and verified to nest the symmetric t exactly at zero skew. A final section turns the tools on data: survival curves for reading tail weight off a log scale, and the Hill estimator recovering the tail index from a sample.
View example →Extreme Value, Survival & Reliability
Eleven laws answering two questions that share a mathematics — how large is the worst case, and when does it fail. Gumbel, Weibull and Fréchet are the three possible shapes for a maximum, and the Generalized Extreme Value law absorbs all three into a single shape parameter ξ; its threshold twin, the Generalized Pareto, governs exceedances. Neither limit theorem is merely asserted — Fisher–Tippett–Gnedenko and Pickands–Balkema–de Haan are both demonstrated by simulation, with block maxima from different parents converging to their GEV and a mean-excess plot going linear. The survival half re-expresses everything through the hazard function h = f/S, where the shape of the failure rate is what classifies a lifetime law: monotone for the Weibull, a non-monotone hump for the log-logistic, and exponential ageing plus constant accident risk for Gompertz–Makeham — with the three components summed explicitly into the classic bathtub curve. The actuarial GB2 closes it as a four-parameter umbrella whose special cases are Burr XII, Dagum and the log-logistic, verified to machine precision rather than by resemblance.
View example →Bayesian Priors & Data Augmentation
The catalog's close, and the one folder about computation rather than data — the eight laws that run inside a modern MCMC sampler. Three jobs organise them: constrain a parameter (the Truncated Normal of probit augmentation, and the Half-Normal / Half-Cauchy / Half-t scale priors), mix scales so a heavy-tailed or sparse prior becomes conditionally Gaussian (the Inverse-Gaussian and the GIG umbrella, whose limits are verified to reproduce the Gamma and Inverse-Gaussian exactly), and linearise a hard likelihood — the job of the Pólya-Gamma, the most consequential modern augmentation, which makes a logistic or negative-binomial likelihood exactly Gaussian given the latent. Having no closed-form density, it is sampled through its infinite Gamma-sum representation and validated through its series density, mean and Laplace transform; the Polson–Scott–Windle identity itself is evaluated on a grid rather than asserted. The Dirichlet closes the set on the simplex. The augmentation demos tie it together — the Bayesian LASSO recovered as a Normal-variance mixture, and the logistic recovered as a Pólya-Gamma mixture.
View example →Multivariate Distributions
The first folder to describe a vector jointly, which introduces the one object the univariate catalog never had — a dependence structure. Six laws: the multivariate Normal and Student-t, elliptical because they depend on x only through the Mahalanobis quadratic form, so their contours are ellipses whose orientation and elongation are the covariance; the skew-normal, which bends them; the Multinomial and its over-dispersed compound the Dirichlet-Multinomial, compared at the same mean so the extra spread is attributable to compounding alone; and the Gaussian Markov random field, where a CAR prior is specified through a sparse precision matrix whose zeros are conditional independences — a dense covariance from a sparse precision, which is why spatial models scale. Every continuous law is drawn twice, as a 3-D density surface and as its iso-density contours, so the dependence is visible rather than merely tabulated. The Gaussian section derives marginals and conditionals in closed form — the two formulas every Gaussian Gibbs update actually is — and a closing copula section strips the margins away to leave dependence on its own.
View example →Matrix-Variate & Covariance Distributions
Where the previous folder treated random vectors, this one treats random matrices — the positive-definite covariance and correlation matrices that are the dependence structure, and the priors you place on them when the covariance itself is unknown. The Wishart is sampled by the Bartlett decomposition, which never forms an outer product; inverting it gives the Inverse-Wishart, the matrix generalisation of the Inverse-Gamma; pairing that with a Gaussian mean gives the Normal-Inverse-Wishart, whose marginal for μ is a heavy-tailed multivariate t — the estimation-risk premium for not knowing Σ. The LKJ prior closes it, putting a single interpretable shape η on correlation matrices, sampled by C-vine and checked against its exact off-diagonal marginal. The hard part of a distribution over matrices is seeing it, so each sampled covariance is drawn as a confidence ellipse (a distribution becomes a cloud of ellipses), correlation matrices become heatmaps, and the space of 2×2 covariances is scattered into its positive-definite cone with every draw confirmed inside the boundary. A closing section builds the modern recommended prior — LKJ correlations times independent Half-Cauchy scales — with the two knobs visibly doing separate work.
View example →Compound & Nonparametric Distributions
The close of the catalog, and the folder where the parameter stops being a fixed vector. A compound law sums a random number of random jumps: the Tweedie adds Gamma jumps a Poisson number of times, producing an atom at zero and a continuous positive part — the exact shape of aggregate insurance claims — with variance following the power law φμp. A nonparametric prior goes further and randomises an entire distribution: stick-breaking snaps a unit stick by repeated Beta(1, α) fractions to give the Dirichlet process, whose clustering view is the Chinese Restaurant Process and whose cluster count grows like α log n — the prior behind mixtures that learn how many components they need. And a Gaussian process randomises an entire function, coupling function values through a kernel so that any finite set of inputs is multivariate Normal, with a closed-form posterior whose uncertainty band tightens at the data and widens in the gaps. None of these has an ordinary density, so each is validated through its construction and moments instead — simulated against exact expectations — and then drawn, since a process is something you look at rather than tabulate.
View example →