A probability distribution belongs to the exponential family if it can be written as p(y;η) = b(y)exp(ηᵀT(y) - a(η)). Here η ∈ ℝᵈ is the natural parameter vector, T(y) ∈ ℝᵈ is the sufficient statistic vector, a(η) is the log-partition function, and b(y) is the base measure. This form is more general than it first appears—by carefully choosing these components, nearly all standard distributions (Normal, Exponential, Chi-square, Poisson, Bernoulli, Beta, Gamma) fit naturally into this framework.

The log-partition function a(η) acts as the normalizing constant. It satisfies ∫b(y)exp(ηᵀT(y))dy = exp(a(η)), ensuring the distribution integrates to one. Computing a(η) from the density reveals its special property: the gradient ∇a(η) = E[T(y)], and the Hessian ∇²a(η) = Cov(T(y)). This means the moments are completely determined by derivatives of a single function—a remarkable compression of information.

Canonical Form and Parameters

The natural parameter η is special: it encodes all the information about the distribution that cannot be separated into the base measure b(y). For a scalar exponential family (d=1), each value of η determines a unique distribution. The mean parameter μ = E[y] can be expressed as μ = da(η)/dη. There is a one-to-one correspondence between natural parameters η and mean parameters μ when a(η) is strictly convex.

The inverse relationship μ⁻¹(η) = da(η)/dη can sometimes be inverted to write η as a function of μ. This relationship is crucial: it tells us how to parameterize the exponential family either by η or by μ. For regression, we typically specify the mean μ = E[y|x] directly and let the model infer η through the exponential family relationship. Understanding this duality between natural and mean parameters is essential for understanding GLMs.

Sufficient Statistics and Data Reduction

The sufficient statistic T(y) summarizes the data: once you know T(y), the actual value of y contains no additional information about the parameters. This allows dramatic data reduction. For Poisson with n observations, the sufficient statistic is the sum Σyᵢ—a single number—yet it contains all information needed to estimate λ. For Bernoulli with n observations, the sufficient statistic is Σyᵢ, the count of successes. For Gaussian with known variance, the sufficient statistic is the sample mean.

This parsimony has profound implications: you do not need to store raw data to do inference—only the sufficient statistics. For streaming data or privacy-sensitive applications, you can compute and transmit sufficient statistics instead of raw observations. Maximum likelihood estimation typically works directly with sufficient statistics, not individual data points. This is why exponential families are so computationally elegant.

Common Examples

Bernoulli (y ∈ {0,1}): p(y;φ) = φʸ(1-φ)¹⁻ʸ gives η = log(φ/(1-φ)), T(y) = y, a(η) = log(1+eⁿ), b(y) = 1. Poisson (y ∈ {0,1,2,...}): p(y;λ) = e⁻λλʸ/y! gives η = log(λ), T(y) = y, a(η) = eⁿ, b(y) = 1/y!. Gaussian (y ∈ ℝ): p(y;μ,σ²) = (1/√(2πσ²))exp(-(y-μ)²/(2σ²)) gives η = μ/σ², T(y) = y, a(η) = η²σ²/2, b(y) = (2πσ²)⁻¹/². Each distribution reveals a unique structure when written in exponential family form.

The Bernoulli distribution models binary random variables y ∈ {0,1}. Its standard form is p(y;φ) = φʸ(1-φ)¹⁻ʸ where φ ∈ (0,1) is the probability of success. To convert this to exponential family form, we algebraically rearrange: p(y;φ) = (1-φ)exp(y log(φ/(1-φ))). This reveals η = log(φ/(1-φ)), the log-odds. The sufficient statistic is T(y) = y, and the log-partition function is a(η) = log(1 + eⁿ).

Notice that a single parameter φ becomes the scalar natural parameter η. The relationship φ = 1/(1 + e⁻ⁿ) is the sigmoid or logistic function. This is not a choice we make when building logistic regression—it emerges automatically from the exponential family structure. The exponential parameterization φ = eⁿ/(1 + eⁿ) normalizes probabilities to [0,1]. The mean is E[y] = φ = da(η)/dη = eⁿ/(1 + eⁿ), confirming the fundamental relationship between natural and mean parameters.

Deriving Logistic Regression

Now assume that the natural parameter varies with features: η = θᵀx. Then E[y|x] = σ(θᵀx) where σ(z) = 1/(1 + e⁻ᶻ). This is logistic regression. The probability of class 1 given x is σ(θᵀx); the probability of class 0 is 1 - σ(θᵀx). The likelihood for a single observation is p(yᵢ|xᵢ,θ) = σ(θᵀxᵢ)^yᵢ (1 - σ(θᵀxᵢ))^(1-yᵢ). The log-likelihood is L(θ) = Σᵢ[yᵢ log σ(θᵀxᵢ) + (1-yᵢ) log(1 - σ(θᵀxᵢ))], which is the cross-entropy loss.

Maximizing this log-likelihood is equivalent to minimizing the cross-entropy. The gradient of the loss is ∇L = Σᵢ(σ(θᵀxᵢ) - yᵢ)xᵢ, which is zero at the maximum likelihood estimate. Unlike linear regression, there is no closed-form solution; we must use gradient descent, Newton's method, or other iterative optimization. The Hessian is ∇²L = -XᵀDX where D is diagonal with entries σ(θᵀxᵢ)(1-σ(θᵀxᵢ)). Since this is negative semidefinite everywhere, the log-likelihood is strictly concave, guaranteeing a unique global maximum.

Interpretation and Odds Ratios

The parameter θⱼ has a natural interpretation via odds ratios. If we increase xⱼ by one unit, holding others fixed, the odds of y=1 change by a factor of exp(θⱼ). If θⱼ = 0.5, increasing xⱼ by one unit multiplies the odds by e^0.5 ≈ 1.65. This makes logistic regression highly interpretable: practitioners can directly translate coefficients into odds ratios without computing probabilities. This interpretability, combined with theoretical grounding in exponential families, makes logistic regression a workhorse for binary classification.

Comparison to Linear Probability Models

An alternative approach would be to assume E[y|x] = θᵀx directly (linear probability model). But this allows predictions outside [0,1], violates the homoscedasticity assumption, and ignores the natural boundary of probability space. The exponential family derivation shows that assuming Bernoulli likelihood and using the canonical link avoids all these issues. The logistic link is not an arbitrary convenience—it is what the mathematics demands.

The multivariate Gaussian distribution is a member of the exponential family. For simplicity, consider the univariate case with fixed variance σ²: p(y;μ) = (1/√(2πσ²))exp(-(y-μ)²/(2σ²)). Expanding the quadratic: -(y-μ)²/(2σ²) = -y²/(2σ²) + μy/σ² - μ²/(2σ²). Rearranging gives p(y;μ) = (1/√(2πσ²))exp(-(y²/(2σ²)) + (μ/σ²)y - μ²/(2σ²)). In exponential family form: η = μ/σ², T(y) = y, a(η) = η²σ²/2, b(y) = (2πσ²)⁻¹/²exp(-y²/(2σ²)).

The natural parameter η equals μ/σ². The log-partition function a(η) = (μ²)/(2σ²) = η²σ²/2. The mean is E[y] = da(η)/dη = η·σ² = μ, recovering the original parameterization. The variance is Var(y) = d²a(η)/dη² = σ², which is constant and not a function of η—unlike Bernoulli where variance depends on the mean. This constant variance assumption is crucial to linear regression.

Deriving Ordinary Least Squares

Assume y|x ~ N(θᵀx, σ²) with σ² fixed. Then E[y|x] = θᵀx, giving the identity link function g(μ) = μ. The likelihood for a single observation is p(yᵢ|xᵢ,θ) = (1/√(2πσ²))exp(-(yᵢ - θᵀxᵢ)²/(2σ²)). The log-likelihood is L(θ) = -n/2 log(2πσ²) - 1/(2σ²)Σᵢ(yᵢ - θᵀxᵢ)². Ignoring constants, minimizing -L(θ) means minimizing Σᵢ(yᵢ - θᵀxᵢ)², which is the sum of squared residuals.

This is ordinary least squares. The maximum likelihood estimate is θ̂ = (XᵀX)⁻¹Xᵀy, derived by setting ∇L = 0. The Hessian is ∇²L = -2XᵀX/σ², which is negative semidefinite when X has full rank, confirming a unique global maximum. The closed-form solution makes OLS computationally efficient and analytically tractable.

Interpreting Coefficients

In linear regression, θⱼ represents the change in expected response per unit increase in xⱼ, holding others constant. If θⱼ = 2.5, then E[y|x] increases by 2.5 for each additional unit of xⱼ. This linear, additive interpretation is transparent and actionable. Confidence intervals for θⱼ can be computed analytically using the standard error SE(θ̂ⱼ) = σ√[(XᵀX)⁻¹ⱼⱼ]. This statistical machinery makes linear regression powerful for inference and prediction.

Assumption Violations

OLS is optimal when the Gaussian assumption holds with constant variance. When y|x has non-Gaussian distribution (e.g., binary, count) or heteroscedastic variance, OLS is no longer the maximum likelihood estimator. Logistic regression is preferred for binary outcomes; Poisson regression for counts. Weighted least squares handles heteroscedasticity. The GLM framework makes clear what to do in each scenario: choose the exponential family matching the data type, then select the canonical link.

A generalized linear model is defined by three key components: (1) a random component specifying the distribution of y|x from the exponential family; (2) a systematic component η = θᵀx, the linear predictor; (3) a link function g relating the mean to the linear predictor via μ = g⁻¹(η). These three components are independent choices—you select the distribution, the link, and the linear model separately, then combine them. This modularity makes GLMs remarkably flexible.

The random component captures how response values scatter around their mean given features. Choosing Bernoulli means y ∈ {0,1}; choosing Gaussian means y ∈ ℝ. The systematic component is always linear in parameters: η = θ₀ + θ₁x₁ + ... + θₚxₚ. The link function transforms the linear predictor into the mean space. The canonical link for an exponential family always satisfies g(μ) = η, making η = θᵀx directly. Using the canonical link simplifies maximum likelihood estimation and ensures the log-likelihood is concave.

Canonical Links

For Bernoulli: μ = φ ∈ (0,1), η = log(φ/(1-φ)), canonical link is logit g(μ) = log(μ/(1-μ)). For Gaussian: μ = E[y] ∈ ℝ, η = μ, canonical link is identity g(μ) = μ. For Poisson: μ = E[y] ∈ (0,∞), η = log(μ), canonical link is log g(μ) = log(μ). For Gamma: μ = E[y] ∈ (0,∞), η = -1/μ, canonical link is reciprocal g(μ) = 1/μ. These emerge from the exponential family form—they are not arbitrary choices. Using canonical links is almost always preferred because they make optimization efficient and interpretable.

Non-Canonical Links

Sometimes domain knowledge or data suggest a non-canonical link. For Bernoulli, the probit link g(μ) = Φ⁻¹(μ) (inverse normal CDF) is sometimes used instead of logit. Both produce predictions in [0,1] and have similar shapes. The main difference is how probabilities near 0.5 relate to the linear predictor: probit has heavier tails. For Poisson, an identity link g(μ) = μ can be used instead of log, though it requires constraints to keep μ > 0. Non-canonical links introduce nonlinearity into optimization but can fit data better if the data violates exponential family assumptions.

Maximum Likelihood Estimation

For any GLM, the log-likelihood is L(θ) = Σᵢ[yᵢθᵢ - a(θᵢ) + log b(yᵢ)] where θᵢ = g(θᵀxᵢ). Taking derivatives: ∂L/∂θⱼ = Σᵢ(yᵢ - μᵢ)xᵢⱼ/Var(yᵢ) where μᵢ = E[yᵢ|xᵢ]. Setting to zero gives the score equation. For canonical links, this simplifies to Σᵢ(yᵢ - μᵢ)xᵢⱼ = 0, which has the intuitive interpretation: weighted residuals are uncorrelated with features. The Hessian is negative definite under mild conditions, ensuring unique global maximum. Numerical methods like iteratively reweighted least squares (IRLS) or gradient descent efficiently find the maximum likelihood estimate.

Advantages of the GLM Framework

The GLM framework unifies diverse regression and classification problems under one mathematical structure. Understanding one GLM helps understand them all. The modular structure (distribution + link + linear model) lets you tailor the model to your data type: choose exponential family matching the response type, choose link matching domain knowledge. The exponential family structure guarantees convex optimization and efficient maximum likelihood estimation. Hypothesis testing and confidence intervals follow standard asymptotic theory. This principled framework beats ad-hoc approaches for credibility, interpretability, and statistical efficiency.

Linear regression emerges from a GLM when we specify: (1) y|x ~ N(μ, σ²) with σ² known and constant; (2) link function g(μ) = μ (identity link); (3) linear predictor η = θᵀx. Since the canonical link for Gaussian is identity, we have η = μ = θᵀx. The predictions are unbounded, which is appropriate for continuous responses. Fitting the model means finding θ that maximizes the likelihood, which is equivalent to minimizing the sum of squared residuals.

The normal equations are XᵀXθ = Xᵀy, solved via θ̂ = (XᵀX)⁻¹Xᵀy when X has full column rank. This closed-form solution is one reason linear regression remains so popular. The predicted values are ŷ = Xθ̂, which lie in the column space of X. Residuals are e = y - ŷ, and they satisfy Xᵀe = 0 (orthogonal to the feature space). The residual sum of squares RSS = eᵀe is minimized by construction. The fitted values form orthogonal projections of y onto the column space of X.

Statistical Properties

Under the Gaussian assumption, the maximum likelihood estimate θ̂ is also the best linear unbiased estimator (BLUE) by the Gauss-Markov theorem. The variance-covariance matrix of θ̂ is Cov(θ̂) = σ²(XᵀX)⁻¹. Unbiased estimation of σ² is σ̂² = RSS/(n - p) where n is sample size and p is number of features. Standard errors SE(θ̂ⱼ) = σ̂√[(XᵀX)⁻¹ⱼⱼ] allow hypothesis tests and confidence intervals. The t-statistic for θⱼ is tⱼ = θ̂ⱼ/SE(θ̂ⱼ), approximately t_{n-p} distributed under the null θⱼ = 0.

Confidence intervals are θ̂ⱼ ± t_{n-p,α/2} · SE(θ̂ⱼ). This analytical inference is a huge advantage of linear regression. Practitioners can quantify uncertainty in predictions and test hypotheses without simulation or bootstrapping. The R² coefficient of determination R² = 1 - RSS/TSS measures goodness of fit, where TSS = Σᵢ(yᵢ - ȳ)² is total sum of squares. A higher R² means more variance explained, though adjustments for model complexity (adjusted R², AIC, BIC) often provide better comparison across models.

When Linear Regression Fails

Linear regression assumes constant variance Var(y|x) = σ² independent of x. When variance changes with x (heteroscedasticity), the estimate θ̂ remains unbiased but standard errors are incorrect, invalidating hypothesis tests. Weighted least squares θ̂ = (XᵀWX)⁻¹XᵀWy with weight matrix W can handle heteroscedasticity. When y is binary, linear regression produces predictions outside [0,1], is less efficient than logistic regression, and violates the homoscedasticity assumption. For count data (non-negative integers), Poisson regression is more appropriate than linear regression. The GLM framework makes clear which model to use: match the exponential family to the response distribution.

Regularization and Shrinkage

When XᵀX is near-singular (multicollinearity), the maximum likelihood estimate becomes unstable and has high variance. Ridge regression adds a penalty: θ̂ᵣᵢdgₑ = (XᵀX + λI)⁻¹Xᵀy. The ridge penalty λ > 0 shrinks coefficients toward zero, trading bias for reduced variance. As λ → 0, ridge approaches OLS; as λ → ∞, all coefficients shrink to zero. Cross-validation selects λ. Lasso adds an ℓ¹ penalty that can drive some coefficients exactly to zero, performing feature selection. Elastic net combines ridge and lasso penalties. These regularization techniques extend linear regression to high dimensions and noisy data, maintaining the GLM framework's interpretability while improving robustness.

Logistic regression is the GLM for binary classification: (1) y|x ~ Bernoulli(φ) with y ∈ {0,1}; (2) link function g(φ) = log(φ/(1-φ)) (logit link); (3) linear predictor η = θᵀx. The predicted probability is φ̂ = σ(θᵀx) = 1/(1 + exp(-θᵀx)). The decision boundary is where φ̂ = 0.5, equivalently where θᵀx = 0. The model's confidence in its prediction increases with |θᵀx|: when |θᵀx| is large, φ̂ is near 0 or 1; when near 0, φ̂ is near 0.5 (maximum uncertainty).

The log-likelihood for a single observation is ℓᵢ(θ) = yᵢ log(σ(θᵀxᵢ)) + (1-yᵢ)log(1-σ(θᵀxᵢ)). Summing over all observations gives the log-likelihood L(θ) = Σᵢℓᵢ(θ). This is the cross-entropy loss, a fundamental objective in machine learning. The gradient is ∇L(θ) = Σᵢ(yᵢ - σ(θᵀxᵢ))xᵢ, which is zero at the MLE. The weighted Hessian is ∇²L(θ) = -XᵀDX where D is diagonal with entries σ(θᵀxᵢ)(1-σ(θᵀxᵢ)). This is negative semidefinite, confirming the log-likelihood is concave and the MLE is unique (up to identifiability issues).

Solving Logistic Regression

Unlike linear regression, logistic regression has no closed-form solution. We use iterative algorithms: Newton's method (most common), gradient descent, or coordinate descent. Newton's method updates θ ← θ - (∇²L)⁻¹∇L, where the inverse Hessian is weighted least squares. Each iteration solves a weighted linear regression problem with weights wᵢ = σ(θᵀxᵢ)(1-σ(θᵀxᵢ)). This is iteratively reweighted least squares (IRLS). Most software packages converge to the MLE in 5-20 iterations. Gradient descent works well with adaptive learning rates like Adam, especially for high-dimensional data.

Computing the standard error of θ̂ⱼ uses the inverse Hessian: SE(θ̂ⱼ) = √[(∇²L)⁻¹ⱼⱼ]. The Wald test statistic zⱼ = θ̂ⱼ/SE(θ̂ⱼ) is approximately N(0,1) under the null θⱼ = 0. Two-sided p-values are 2Φ(-|zⱼ|). Confidence intervals are θ̂ⱼ ± 1.96·SE(θ̂ⱼ) at the 95% level. These asymptotic results are accurate for large n, but may be poor for small samples, especially when outcomes are rare. Penalized maximum likelihood (ridge, lasso) improves stability, especially when features are highly correlated.

Multiclass Extension

Binary logistic regression extends to K > 2 classes via softmax regression. Let πₖ = P(y=k|x) for k=0,...,K-1. Then log(πₖ/π₀) = θₖᵀx for k=1,...,K-1. The model learns K-1 parameter vectors θ₁,...,θₖ₋₁ (setting θ₀=0 for identifiability). The predicted probabilities are πₖ = exp(θₖᵀx)/Σⱼexp(θⱼᵀx), where we set θ₀=0. This is the softmax function. One-vs-rest and one-vs-one approaches are alternatives, but they lack the principled probabilistic interpretation of multinomial logistic regression.

Interpretation and Odds Ratios

If θⱼ = 0.5, increasing xⱼ by one unit (holding others fixed) multiplies the odds of y=1 by exp(0.5) ≈ 1.65. Conversely, if θⱼ = -0.2, increasing xⱼ by one unit divides the odds by exp(0.2) ≈ 1.22, reducing the odds of y=1. This odds ratio interpretation makes logistic regression highly actionable: coefficients translate directly into business impact without computing probabilities. The linearity in log-odds is often reasonable for many domains, though interactions and nonlinearities can be added via polynomial or spline features.