Key Distributions for ML

A product manager asks why your demand-forecasting model outputs a point estimate plus an interval, and why the interval is symmetric. The honest answer traces to one modeling decision: you assumed the residuals are normally distributed. The normal (Gaussian) distribution is not the default in ML by tradition alone — it is the unique distribution that shows up wherever many small independent effects add together, which is why it appears in regression noise, weight initialization, and the central limit theorem.

The normal is parameterized by mean μ (location) and standard deviation σ (scale): X ~ N(μ, σ²). Its density is the familiar bell curve, f(x) = (1/(σ√(2π)))·exp(−(x−μ)²/(2σ²)). Two properties carry most of its ML mileage:

1. The 68–95–99.7 rule. ~68% of mass lies within μ ± σ, ~95% within μ ± 2σ, ~99.7% within μ ± 3σ. This is what makes a symmetric prediction interval computable without simulation: μ ± 1.96σ covers 95% of a normal's mass.
2. Stability under sums and scaling. A sum of independent normals is normal with means and variances adding: N(μ₁, σ₁²) + N(μ₂, σ₂²) = N(μ₁+μ₂, σ₁²+σ₂²). Standardizing any normal via Z = (X − μ)/σ gives the single canonical N(0, 1), which is why z-scores and standard scaling exist at all.

Where ML leans on it:

The likelihood connection deserves one line of algebra, because it is the deepest item on that table. If y = f(x) + ε with ε ~ N(0, σ²), then maximizing the log-likelihood of the data is minimizing Σ(yᵢ − ŷᵢ)² — ordinary least squares is literally negative-Gaussian-log-likelihood. Switch the noise assumption to another family and the loss changes with it; this is why L1/Laplace noise yields MAE instead of MSE.

from scipy.stats import norm
import numpy as np

# Forecast errors: mu = 0, sigma = 120 units/day
mu, sigma = 0, 120
print(norm.cdf(1.96 * sigma, loc=mu, scale=sigma) - norm.cdf(-1.96 * sigma, loc=mu, scale=sigma))
# 0.95 — the 95% interval

# Xavier-style init for a layer with 256 fan-in
rng = np.random.default_rng(0)
w = rng.normal(0, np.sqrt(1 / 256), size=(256, 256))
print(w.mean(), w.std())  # ≈ 0, ≈ 0.0625 — activations keep their scale

An interviewer's favorite trap: "the data is normal, right?" Real features like revenue, latency, and counts are right-skewed — the normal is a modeling assumption you check and often replace (log transform), not a law. You measured exactly that with skewness in Descriptive Statistics and EDA; the Bernoulli and binomial families next handle the binary targets where the normal plainly does not apply.

Your ad model predicts, per impression, whether a user clicks. One impression has exactly two outcomes — click (1) or no click (0) — with some underlying probability p. That single trial is a Bernoulli(p) variable: P(X = 1) = p, P(X = 0) = 1 − p, mean p, variance p(1−p). Note the variance is not free: the noisier the coin (p near 0.5), the higher the variance, peaking at 0.25 and vanishing as p → 0 or 1.

Repeat the trial n times independently and count successes: X ~ Binomial(n, p), with P(X = k) = C(n,k)·pᵏ·(1−p)ⁿ⁻ᵏ. Its parameters are the natural summary of repeated binary events:

This is the moment to see the whole classification stack as one machine. Logistic regression models log-odds of a Bernoulli probability: it learns w so that σ(wᵀx) ≈ p = P(click | features), where σ is the sigmoid. The binary cross-entropy loss is exactly the negative log-likelihood of Bernoulli(p = σ(wᵀx)) observed at label y ∈ {0, 1}: −[y·log p + (1−y)·log(1−p)]. When someone asks "why not square error on classification targets?", this is the answer: square error is the Gaussian-noise loss from the previous chapter; labels here are Bernoulli draws, so the matching likelihood is Bernoulli/cross-entropy.

Two modeling caveats interviewers probe. First, the binomial assumes independence with a fixed p across trials. Real ad impressions violate this: the same user sees three ads in one session (positively correlated, so the effective sample size is smaller than n) and per-user click rates differ (heterogeneous p). Both violations inflate the variance beyond np(1−p) — a condition called overdispersion — and the standard fix is a beta-binomial or a mixed-effects model. Second, np is fine to interpret only when n is large and p not tiny; otherwise use the exact PMF, not the normal approximation.

from scipy.stats import binom, bernoulli

n, p = 10, 0.30
print(binom.mean(n, p), binom.var(n, p))        # 3.0, 2.1
print(binom.pmf(3, n, p))                       # 0.2668 — P(exactly 3 clicks)
print(binom.pmf(range(n + 1), n, p).round(3))   # full PMF, sums to 1

# Sanity-check a week of logged data against the assumed binomial
observed_clicks = [2, 4, 3, 1, 5, 3, 2]        # per-impression-batch counts, n=10
mean_obs, var_obs = np.mean(observed_clicks), np.var(observed_clicks, ddof=1)
print(mean_obs, var_obs)
# var_obs >> np(1-p)=2.1 would flag overdispersion → revisit independence

# Logistic regression is Bernoulli likelihood in disguise
from sklearn.linear_model import LogisticRegression
m = LogisticRegression().fit(X_train, y_train)  # minimizes Bernoulli neg-log-likelihood

The empirical overdispersion check above — compare the observed variance of counts to np(1−p) — is the fastest way to honor the distribution's assumptions instead of asserting them. When counts arrive over time or space rather than as batched trials, the right family changes again, and that is the Poisson's job.

You are sizing servers for an API and have a table of request counts per minute. The values are non-negative integers, there is no fixed number of "trials," and the events arrive continuously over time. The normal is a poor fit (it assigns mass to negative counts and treats the data as continuous); the binomial is unavailable (no n). The Poisson(λ) distribution is purpose-built for this: X = number of events in a fixed interval, with P(X = k) = e⁻λ·λᵏ/k!, mean λ, and variance λ.

The signature fact — mean equals variance — is the fastest diagnostic you will ever run on count data. If the empirical variance is roughly double the empirical mean (variance-to-mean ratio ≈ 2), the Poisson assumption is already dead: real count processes (requests, defects, tweets) usually exhibit overdispersion because "λ itself fluctuates" with traffic bursts. The fix is the negative binomial, the count-data analogue of what you saw with binomial overdispersion in the previous chapter.

Where each family earns its place:

A small λ is what makes Poisson a rare event model: with λ = 0.5 requests/minute, P(X ≥ 3) = 1 − P(0) − P(1) − P(2) ≈ 0.014 — that small tail probability is exactly what capacity planning needs to know. And under a Poisson process, counts in disjoint time intervals are independent, which reconnects to the independence machinery from Probability Foundations.

from scipy.stats import poisson
import numpy as np

lam = 0.5
print(poisson.pmf([0, 1, 2, 3], lam).round(4))   # [0.6065 0.3033 0.1008 0.0168]
print(1 - poisson.cdf(2, lam))                    # 0.0144 — P(≥3 in one minute)

# Poisson regression: log of the mean is linear in the features
import statsmodels.api as sm
model = sm.GLM(counts, features,
               family=sm.families.Poisson()).fit()  # models log(λ) = wᵀx

# Mean-vs-variance diagnostic on real logs
lam_hat = counts.mean()
print(counts.mean(), counts.var(ddof=1), counts.var(ddof=1) / lam_hat)
# ratio far above 1 → overdispersion → switch to negative binomial

The GLM call shows the modern role of the Poisson in ML: Poisson regression constrains predictions to be non-negative counts and models log(λ) linearly, which beats a plain regressor that can predict −3 server requests. Note also how λ enters only through log λ in the model — the distribution choice, not just the data type, decided the link function.

You now have the three workhorse families and the discipline of checking their assumptions against empirical moments. The next topic turns that discipline into a procedure: Bayesian Thinking for Models uses these distributions as likelihoods and updates beliefs about their parameters as data arrives.