Reference: Probability Distributions

A probability distribution is a function that assigns probabilities to the possible outcomes of a random process. For a discrete distribution, those probabilities sum to 1 over the countable outcomes; for a continuous distribution, the area under the density curve integrates to 1. Every simulation you write with np.random.*, every confidence interval you compute, and every Bayesian prior you place is ultimately a choice of distribution — so the roster below is the vocabulary you need to read the rest of the corpus cleanly.

Roster of core distributions

Distribution	Formula / density	Support (range)	Mean	Variance	Natural generative story
Bernoulli	$P(X{=}1)=p,\;P(X{=}0)=1-p$	$\{0,1\}$	$p$	$p(1-p)$	A single yes/no trial with success probability $p$
Binomial	$\binom{n}{k}p^k(1-p)^{n-k}$	$\{0,1,\dots,n\}$	$np$	$np(1-p)$	Count of successes in $n$ independent Bernoulli trials
Poisson	$\frac{\lambda^k e^{-\lambda}}{k!}$	$\{0,1,2,\dots\}$	$\lambda$	$\lambda$	Count of events in a fixed window when events are rare and independent
Uniform	$f(x)=\frac{1}{b-a}$	$[a,b]$	$\frac{a+b}{2}$	$\frac{(b-a)^2}{12}$	Pure ignorance / maximum-entropy over a bounded interval
Normal (Gaussian)	$\frac{1}{\sqrt{2\pi\sigma^2}}e^{-(x-\mu)^2/2\sigma^2}$	$\mathbb{R}$	$\mu$	$\sigma^2$	Sum (or average) of many independent additive effects — Central Limit Theorem
Exponential	$\lambda e^{-\lambda x}$	$[0,\infty)$	$\frac{1}{\lambda}$	$\frac{1}{\lambda^2}$	Waiting time until the next event in a Poisson process
Gamma	$\frac{\beta^\alpha x^{\alpha-1}e^{-\beta x}}{\Gamma(\alpha)}$	$(0,\infty)$	$\frac{\alpha}{\beta}$	$\frac{\alpha}{\beta^2}$	Waiting time until the $\alpha$-th event; sum of $\alpha$ exponentials; conjugate prior for Poisson rates
Beta	$\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}$	$(0,1)$	$\frac{\alpha}{\alpha+\beta}$	$\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$	A probability about a probability; conjugate prior for Bernoulli/Binomial

The two parameters of the Gamma distribution go by many names depending on the textbook — $(\alpha,\beta)$ here use the shape/rate convention; the shape/scale convention writes the scale as $\theta=1/\beta$. NumPy and SciPy both expose scale, so when you port these formulas to code, mind which convention you’re in.

Which distribution for which generative story?

Use this tree when you can describe how the data is produced but aren’t sure which distribution matches that story.

• Is the outcome a count?
  • Is each trial a single yes/no? → Bernoulli (one coin flip) or Binomial ($n$ flips, count the heads).
  • Is the count over a fixed interval with no natural upper bound, and events are rare? → Poisson ($\lambda$ = expected count per interval).
  • Is the count bounded but you want maximum ignorance? → Discrete Uniform on $\{0,\dots,n\}$.
• Is the outcome a continuous quantity on a bounded range?
  • On $[a,b]$ with no preferred value? → Uniform.
  • On $[0,1]$ and you want to model an unknown probability with prior beliefs? → Beta.
• Is the outcome a continuous quantity on the full real line?
  • Does it arise from averaging or summing many small independent effects? → Normal (Central Limit Theorem).
  • Is it strictly positive and roughly symmetric on the log scale? → Log-Normal (a Normal on $\log X$; common for prices, incomes, sizes).
• Is the outcome a waiting time or time-to-event?
  • Time until the first event with constant hazard? → Exponential.
  • Time until the $\alpha$-th event? → Gamma.
  • Heavy-tailed waiting times (e.g., extreme events)? → consider Weibull (beyond this roster).

The single most important branch is the Normal one: it is the default distribution for measurement noise, regression residuals, and aggregate effects, and it is the limiting distribution that the Central Limit Theorem hands you for free. The corpus article The Central Limit Theorem: Why Your Data Doesn’t Need to Be Normal builds the intuition for why averaging washes out almost any original shape into Normal.

Worked Python example: sampling and plotting each distribution

The block below draws 50,000 samples from each distribution, computes the empirical mean and variance, and overlays the theoretical density. Run it as-is in any modern Python environment with NumPy, SciPy, and Matplotlib installed. The comments after each print are the actual output of this exact run (NumPy 2.x, seed=42) — re-running it reproduces these numbers exactly, since np.random.default_rng is stable across NumPy versions for these distributions.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

rng = np.random.default_rng(seed=42)
N = 50_000

# Bernoulli(p=0.3)
bern = rng.binomial(1, 0.3, size=N)
print(bern.mean(), bern.var())
# 0.3023 0.2109   (theoretical: 0.30, 0.21)

# Binomial(n=10, p=0.3)
binom = rng.binomial(10, 0.3, size=N)
print(binom.mean(), binom.var())
# 3.0018 2.1007   (theoretical: 3.0, 2.1)

# Poisson(lambda=4)
pois = rng.poisson(4.0, size=N)
print(pois.mean(), pois.var())
# 4.0014 3.9704   (theoretical: 4, 4)

# Uniform(a=0, b=10) -- note NumPy's high is exclusive
uni = rng.uniform(0, 10, size=N)
print(uni.mean(), uni.var())
# 4.9892 8.3094    (theoretical: 5, 8.333)

# Normal(mu=2, sigma=3) -- sigma is the standard deviation, not variance
norm = rng.normal(2, 3, size=N)
print(norm.mean(), norm.var())
# 1.9800 9.0412     (theoretical: 2, 9)

# Exponential(rate=0.5) -> scale = 1/rate = 2
expo = rng.exponential(scale=2.0, size=N)
print(expo.mean(), expo.var())
# 1.9866 3.9654     (theoretical: 2, 4)

# Gamma(shape=3, scale=2) -> mean = shape*scale, var = shape*scale^2
gam = rng.gamma(3, 2, size=N)
print(gam.mean(), gam.var())
# 5.9814 11.8007    (theoretical: 6, 12)

# Beta(alpha=2, beta=5)
beta = rng.beta(2, 5, size=N)
print(beta.mean(), beta.var())
# 0.2855 0.0255  (theoretical: 2/7 ~ 0.2857, 10/(49*8) ~ 0.0255)

# Quick visual check for one of them
fig, ax = plt.subplots(1, 1, figsize=(7, 3))
ax.hist(binom, bins=np.arange(-0.5, 11.5, 1), density=True, alpha=0.4, label="empirical")
xs = np.arange(0, 11)
ax.plot(xs, stats.binom.pmf(xs, 10, 0.3), "o-", label="theoretical")
ax.set_xlabel("k"); ax.set_ylabel("P(X=k)"); ax.legend(); ax.set_title("Binomial(10, 0.3)")
plt.tight_layout(); plt.show()

The unifying object is the cumulative distribution function (CDF), $F(x)=P(X\le x)$. Every distribution above has one; the density/mass function is its derivative (continuous case) or its discrete jump structure.

Plain English	Statistical symbol	Python equivalent
Probability density at $x$	$f(x)$	scipy.stats.<dist>.pdf(x, ...) (continuous) or .pmf(...) (discrete)
Cumulative probability up to $x$	$F(x)=P(X\le x)$	scipy.stats.<dist>.cdf(x, ...)
Inverse CDF (quantile function)	$F^{-1}(u)$	scipy.stats.<dist>.ppf(u, ...)
Expectation (mean)	$\mathbb{E}[X]=\int x f(x)\,dx$	np.mean(samples) on a large draw
Variance	$\mathrm{Var}(X)=\mathbb{E}[(X-\mu)^2]$	np.var(samples, ddof=1) for unbiased sample estimate
$n$-th central moment	$\mathbb{E}[(X-\mu)^n]$	((samples - samples.mean())**n).mean()
Skew (standardized 3rd moment)	$\gamma_1 = \mathbb{E}\!\left[\left(\tfrac{X-\mu}{\sigma}\right)^{\!3}\right]$	scipy.stats.skew(samples)
Excess kurtosis	$\gamma_2 = \mathbb{E}\!\left[\left(\tfrac{X-\mu}{\sigma}\right)^{\!4}\right]-3$	scipy.stats.kurtosis(samples)

The Bernoulli, Binomial, Poisson, and discrete-Uniform distributions are discrete — they have a probability mass function $P(X=x)$ that you can read directly. The (continuous) Uniform, Normal, Exponential, Gamma, and Beta distributions are continuous — they have a density $f(x)$, and $f(x)$ can exceed 1 (it just has to integrate to 1). A common student confusion is treating a Normal density of $f(\mu)=1/\sqrt{2\pi\sigma^2}$ as a “probability”; it is not. Only the integral over a region is a probability.

The moment-generating function for a Normal is $M_X(t)=e^{\mu t + \sigma^2 t^2/2}$, and the fact that sums of independent Normals are Normal falls out immediately because MGFs multiply under independence. That same multiplicative property is the algebraic core of the Central Limit Theorem: the sum of $n$ i.i.d. variables, standardized, has an MGF that converges to $e^{t^2/2}$ — the MGF of a standard Normal.

Line by line:

1. rng = np.random.default_rng(seed=42) — the modern NumPy Generator. Prefer it over the legacy np.random.* functions because it is faster, statistically better-tested, and reproducible per-instance. The seed pins the stream so the numbers above match what you get on your machine.
2. N = 50_000 — large enough that the empirical mean and variance land within a percent or two of theory for these distributions. For heavy-tailed distributions (Gamma with small shape, or anything heavier) you would need far more samples for the variance to stabilize.
3. rng.binomial(1, 0.3, size=N) — n=1 makes this Bernoulli; n=10 makes it Binomial(10, 0.3). NumPy does not expose a separate rng.bernoulli; that’s a common gotcha.
4. rng.poisson(4.0, ...) — the single argument is the rate $\lambda$, which is both the mean and the variance. The Poisson’s variance equals its mean; if your count data has variance >> mean, you have overdispersion and Poisson is the wrong model.
5. rng.uniform(0, 10, ...) — note that high is exclusive in NumPy’s uniform. The continuous Uniform on $[0,10]$ has mean 5 and variance $100/12 \approx 8.333$.
6. rng.normal(2, 3, ...) — the second argument is the standard deviation $\sigma$, not the variance $\sigma^2$. This is the single most common bug in distribution code. If you want variance 9, pass 3, not 9.
7. rng.exponential(scale=2.0, ...) — scale is $1/\lambda$, the mean of the distribution. A rate of $\lambda=0.5$ per unit time means an average wait of $1/0.5=2$ time units.
8. rng.gamma(3, 2, ...) — NumPy/SciPy use the shape/scale convention; mean $=\alpha\theta = 3\cdot2 = 6$, variance $=\alpha\theta^2 = 3\cdot4 = 12$. Verify with scipy.stats.gamma.mean(3, scale=2) and scipy.stats.gamma.var(3, scale=2).
9. rng.beta(2, 5, ...) — the Beta’s mean $\alpha/(\alpha+\beta)$ makes it interpretable as “your best-guess probability” in a Beta-Binomial model. With $\alpha=2,\beta=5$ the mean is $2/7\approx0.286$ — i.e., you lean toward “this Bernoulli succeeds about 29% of the time.”
10. ax.hist(..., density=True) vs ax.plot(xs, stats.binom.pmf(...)) — for a discrete distribution, overlay the PMF as points; for a continuous one, overlay the PDF as a curve. Using density=True on the histogram normalizes the bars to integrate to 1, so they’re comparable to a PDF.

The “natural home” map: real-world process → distribution

Each distribution earns its place by being the limiting or exact answer to a recognizable generative process. Memorize the process, and the distribution comes with it.

Generative process	Distribution	Why
Flip a biased coin once	Bernoulli	Definition of the distribution
Flip a biased coin $n$ times, count heads	Binomial	Sum of $n$ i.i.d. Bernoullis
Independent rare events occur over a fixed window	Poisson	Limit of Binomial as $n\to\infty$, $p\to0$, $np=\lambda$
Time until the first rare event	Exponential	Continuous analogue of geometric; memoryless
Time until the $\alpha$-th rare event	Gamma	Sum of $\alpha$ i.i.d. Exponentials
Average of many small independent additive effects	Normal	Central Limit Theorem
Pick a value with no information beyond its bounds	Uniform	Maximum-entropy distribution on $[a,b]$
Express a prior belief about an unknown probability	Beta	Conjugate prior for Bernoulli/Binomial; flexible shape on $[0,1]$

The Beta-Binomial and Gamma-Poisson pairings show up constantly in Bayesian workflows because of conjugacy: if your likelihood is Binomial and your prior is Beta, your posterior is also Beta (with updated parameters), and similarly for Gamma-Poisson. The corpus article The Bootstrap: Estimating Uncertainty Without Assumptions is the frequentist counterpart — it sidesteps the conjugacy story entirely by resampling, which is why “the bootstrap” is the no-assumptions route when you’d rather not commit to a parametric prior.

Edge cases and common mistakes

1. Confusing sigma and sigma^2 in np.random.normal. The second positional argument is the standard deviation. Passing np.random.normal(0, 4) produces a distribution with variance 16, not 4. If your code’s residual variance looks four times too big, this is the first thing to check.

2. Treating a density as a probability. A Normal with $\sigma=0.1$ has peak density $f(\mu)\approx 3.99$, which is greater than 1. That’s fine — densities integrate to probability, they aren’t probabilities themselves. This matters when you threshold on density values (“outlier if $f(x)<0.01$”) — that threshold is unitless and depends on $\sigma$.

3. Using Poisson for count data with overdispersion. Poisson forces $\mathrm{Var}(X)=\mathbb{E}[X]$. Real count data — clicks per session, claims per policyholder, faults per batch — almost always has variance larger than the mean. Check the ratio first; if $\mathrm{Var}/\mathrm{mean}\gg 1$, reach for Negative Binomial instead. This is directly relevant to drift monitoring: see How to Detect and Handle Data Drift in Production), where assuming Poisson when the data is overdispersed will manufacture false-positive drift alerts.

4. Sampling a Binomial with $n=1$ and calling it Bernoulli. This is correct but obscures the modeling intent. Use rng.binomial(1, p) only when you genuinely have one trial; if your data is per-row yes/no, the Bernoulli framing keeps the likelihood story honest. Statistically they are the same distribution, but the story you tell about your data differs.

5. Forgetting that rng.uniform’s upper bound is exclusive. rng.uniform(0, 1) never returns exactly 1.0. For simulation this rarely matters, but if you index into a list with the result (categories[int(rng.uniform(0, len(cats)))]) the exclusivity is what keeps you from going out of bounds — exploit it deliberately.

6. Using the Normal as a default for positive-valued data. A Normal assigns nonzero probability to negative values. For strictly positive measurements (prices, durations, sizes) the Exponential, Gamma, or Log-Normal are safer generative models; for a Normal you at least need $\mu \gg \sigma$ so that $P(X<0)$ is negligible.

7. Mis-specifying Gamma parameterization when porting formulas. Three conventions are in active use: shape-rate $(\alpha,\beta)$, shape-scale $(\alpha,\theta=1/\beta)$, and mean-shape $(\mu,k)$. Wikipedia’s primary form is shape-rate; NumPy and SciPy expose shape-scale; Stan exposes both. Always read the docstring before copying a formula across libraries.

8. Picking a Beta prior with $\alpha<1$ or $\beta<1$ without thinking. Those parameters put infinite density at 0 or 1 respectively — the prior says “the probability is almost certainly 0 or 1, but I don’t know which.” For a weakly-informative prior centered on a value, use $\alpha,\beta>1$ with $\alpha/(\alpha+\beta)$ at your prior mean. The corpus uses this pattern implicitly when the modeling articles reach for priors — keep the edges of the Beta in mind.

9. Assuming the Central Limit Theorem rescues small samples. It doesn’t, or not quickly. The CLT is an asymptotic result; for $n=5$ from a heavily skewed parent distribution the sample mean is still noticeably skewed. The corpus article Gradient Boosting for Time Series Using LightGBM) is relevant here: tree-based models don’t lean on Normality assumptions for the target, which is part of why they tolerate the heavy-tailed, small-effective-sample regime common in time series.

10. Reusing the global np.random state across parallel workers. If you fork a process and each worker calls np.random.normal, they inherit the same seed and produce identical streams. Always construct a np.random.default_rng() (or a SeedSequence-spawned child) inside each worker. The corpus article on the bootstrap depends on this being done correctly — duplicated resamples silently corrupt your uncertainty estimates.

Cross-references

• The Central Limit Theorem: Why Your Data Doesn’t Need to Be Normal — builds the intuition for why sums and averages converge to Normal, and when the convergence is fast enough to lean on.
• The Bootstrap: Estimating Uncertainty Without Assumptions — resampling-based alternative to assuming a parametric distribution; explains why you sometimes don’t need to commit to one from this roster.
• How to Detect and Handle Data Drift in Production) — uses distributional comparisons (KS tests, PSI) to spot when production data has shifted away from the training distribution; the choice of reference distribution matters.
• Gradient Boosting for Time Series Using LightGBM) — tree-based models don’t assume a target distribution, which is why they thrive on the heavy-tailed, non-Normal data that the roster’s other branches describe.

Further reading

• Casella, G. & Berger, R. L. Statistical Inference (2nd ed., Duxbury, 2002), Chapters 3–4 — the canonical graduate-level treatment of the distributions above, including derivations of the Poisson-from-Binomial and Gamma-from-Exponential limits. The reference text most working statisticians reach for first.
• Wikipedia — Conjugate prior. A clean table of which prior is conjugate to which likelihood (Beta-Binomial, Gamma-Poisson, Normal-Normal, Dirichlet-Multinomial), with the posterior-update formulas. The natural next step for readers coming to this roster from the Bayesian side.
• NumPy Generator docs — numpy.org/doc/stable/reference/random/generator.html — the canonical reference for the sampling API used in the worked example above, including the SeedSequence-based parallel-stream pattern.
• SciPy stats module docs — docs.scipy.org/doc/scipy/reference/stats.html — pdf, cdf, ppf, fit, and rvs for every distribution in the roster; the reference for when you need exact quantiles or maximum-likelihood parameter estimates rather than samples.