Exponential Family & Link Functions — Exponential family and link functions - Statistics Graduate Level

SECTION 1

Historical Context & Motivation

Statistical modeling in the early twentieth century was largely confined to the normal (Gaussian) distribution: ordinary least squares regression assumed that the response variable was continuous and that errors were normally distributed with constant variance. While this framework proved remarkably powerful for many applications, it offered no principled way to handle count data, binary outcomes, or strictly positive continuous responses. Researchers in diverse fields—epidemiology, ecology, actuarial science—needed models for proportions, rates, and durations, yet the tools available forced them into awkward ad hoc transformations or separate theories for each distribution. The quest for a unified regression framework capable of accommodating multiple data types motivated the development of the exponential family formalism and its associated link functions.

1934–36 • Exponential Family Formalized

E. J. G. Pitman and B. O. Koopman independently identified a broad class of distributions sharing a common algebraic form, establishing the concept of sufficient statistics and laying the mathematical groundwork for what would become the exponential family.

1958 • Fisher's Information & MLE Theory

Building on R. A. Fisher's earlier work, the complete theory of maximum likelihood estimation for exponential families was crystallized, proving that the MLE achieves the Cramér–Rao lower bound asymptotically—a property uniquely elegant within this family.

1972 • Nelder & Wedderburn Introduce GLMs

John Nelder and Robert Wedderburn published their landmark paper in JRSS-A, defining generalized linear models (GLMs). By pairing an exponential-family distribution with a link function and a linear predictor, they unified logistic regression, Poisson regression, and classical linear regression into a single algorithmic framework.

1983 • McCullagh & Nelder's GLM Textbook

The publication of 'Generalized Linear Models' by McCullagh and Nelder provided a comprehensive treatment of the theory, establishing the canonical link function concept, deviance-based diagnostics, and the iteratively reweighted least squares (IRLS) algorithm as standard tools.

2000s–present • Modern Extensions

The exponential-family GLM framework now underpins generalized additive models (GAMs), mixed-effects models (GLMMs), Bayesian hierarchical models, and regularized regression (e.g., elastic net for GLMs), serving as the foundational grammar of modern applied statistics.

The central question these developments addressed can be stated concisely: given a response variable whose distribution is not Gaussian—perhaps Bernoulli, Poisson, gamma, or binomial—how can we build a systematic regression model that respects the distributional constraints while still permitting efficient parameter estimation? The answer lies in recognizing that all of these distributions belong to a single algebraic family and that a carefully chosen link function bridges the distribution's natural parameter space with the unconstrained real line on which linear predictors operate.

SECTION 2

Core Principles & Definitions

The theoretical architecture of generalized linear models rests on a small number of interlocking ideas. A probability distribution belongs to the exponential family if its density (or mass) function can be written in a particular canonical form. Once that structure is recognized, the GLM framework specifies three components: a random component (the exponential-family distribution), a systematic component (the linear predictor η = Xβ), and a link function that maps the mean of the response to the linear predictor. Understanding these pieces and how they connect is the core of the subject.

Exponential-Family Form

A distribution belongs to the exponential family if its density can be written as f(y|θ, φ) = a(y, φ) · exp{[yθ − b(θ)] / φ}, where θ is the natural (canonical) parameter, φ is the dispersion parameter, and b(θ) is the log-partition function.

Mean–Variance Relationship

For exponential-family members, the mean is μ = b′(θ) and the variance is Var(Y) = φ · b″(θ). This yields a deterministic mean–variance relationship that distinguishes each distribution—for example, Var(Y) = μ for Poisson and Var(Y) = μ(1 − μ) for Bernoulli.

Link Function g(·)

The link function g maps the mean μ to the linear predictor: g(μ) = η = x⊤β. It must be monotone and differentiable. The canonical link is the special case where g(μ) = θ, producing mathematical simplifications in estimation and inference.

Canonical Link & Sufficiency

When the canonical link is used, X⊤Y is a sufficient statistic for β. This simplifies the score equations to X⊤(Y − μ) = 0 and guarantees concavity of the log-likelihood, ensuring a unique maximum likelihood estimate.

Iteratively Reweighted Least Squares

GLM parameters are estimated via IRLS, which at each iteration solves a weighted least squares problem using a working response and a weight matrix derived from the current mean and variance estimates—effectively a Fisher scoring algorithm.

✦

KEY TAKEAWAY

Think of the exponential family as a universal adapter standard—like USB-C for electronics. Just as USB-C lets a single port connect to monitors, storage devices, and chargers by using different protocols over the same physical interface, the exponential-family form lets a single estimation algorithm (IRLS) handle Bernoulli, Poisson, Gaussian, gamma, and other distributions by routing through different link functions and variance functions. The canonical form is the shared specification that makes everything plug in.

SECTION 3

Visual Explanation — The GLM Architecture

The diagram illustrates the three-component architecture of a GLM. The random component specifies the exponential-family distribution of Y. The link function maps the mean μ onto the real-valued linear predictor η. The systematic component is the familiar linear combination of covariates. The table at the bottom shows how different distributions pair with their canonical links.

The diagram above makes a critical structural point: in a GLM, the modeler's choice reduces to selecting which distribution the data follow (the left box) and which link function connects the mean to the predictor (the center box). The linear predictor on the right remains structurally identical across all GLMs—it is always a linear combination of covariates. This separation of concerns is what makes the framework so powerful. When we use the canonical link, the connection between the natural parameter θ and the linear predictor η becomes an identity, which yields especially clean score equations and guarantees the existence and uniqueness of the maximum likelihood estimator under mild regularity conditions.

SECTION 4

Mathematical Framework

We now develop the exponential-family representation rigorously and derive the key relationships that underpin GLM estimation. All results assume a single-parameter exponential family with a possible dispersion parameter φ, which is the standard setting for GLMs.

EXPONENTIAL-FAMILY DENSITY

f(y | θ, φ) = a(y, φ) · exp{ [y·θ − b(θ)] / φ }

Here θ is the natural (canonical) parameter, φ is the dispersion parameter (φ = 1 for Bernoulli and Poisson), b(θ) is the cumulant (log-partition) function, and a(y, φ) is a normalizing factor that does not depend on θ.

MEAN AND VARIANCE FROM b(θ)

E[Y] = μ = b′(θ) Var(Y) = φ · b″(θ)

These follow from differentiating the log-likelihood. The first derivative of b with respect to θ gives the mean; the second derivative gives the variance function V(μ) up to the dispersion factor. This is the source of the mean–variance relationship that characterizes each member of the family.

To see why these identities hold, consider the log-likelihood for a single observation: ℓ(θ) = log a(y, φ) + [yθ − b(θ)] / φ. Taking the first derivative with respect to θ and setting it to zero at the true parameter gives E[∂ℓ/∂θ] = 0, which implies E[Y] = b′(θ). The second derivative identity Var(∂ℓ/∂θ) = −E[∂²ℓ/∂θ²] then yields Var(Y) = φ · b″(θ). Because b″(θ) > 0 for a well-defined exponential family, the log-likelihood is strictly concave in θ, guaranteeing a unique MLE.

LINK FUNCTION AND LINEAR PREDICTOR

g(μ) = η = x⊤β μ = g⁻¹(η)

The function g is a monotone, differentiable link function. When g(μ) = θ, we call it the canonical link. Non-canonical links (e.g., probit for binary data) are also permissible but sacrifice some algebraic simplicity.

SCORE EQUATION (CANONICAL LINK)

X⊤(Y − μ̂) = 0

Under the canonical link, the score equations reduce to this elegant form, which mirrors the normal equations of ordinary least squares. X⊤Y is a sufficient statistic for β. The system is solved iteratively via IRLS: at iteration t, compute the working response z and weight matrix W from current estimates, then update β by solving (X⊤WX)β = X⊤Wz.

ℹ

Why the Canonical Link is Special

When g(μ) = θ, the natural parameter equals the linear predictor, and the log-likelihood becomes linear in β through the term yθ = y·x⊤β. This linearity ensures that the Fisher information matrix equals the negative Hessian (observed = expected information), the log-likelihood is globally concave in β, and X⊤Y is minimal sufficient. These properties simplify both theory and computation. For non-canonical links, the log-likelihood may still be concave in practice, but this is not guaranteed in general, and the score equations require an additional weight correction involving dμ/dη.

SECTION 5

Exponential-Family Members in Detail

Each common GLM response distribution can be cast into the exponential-family form by identifying θ, b(θ), φ, and a(y, φ). The following diagram provides a visual taxonomy, and the subsequent table offers a comprehensive reference. Recognizing these mappings is an essential skill: when faced with a new distribution, systematically extracting the canonical parameter and cumulant function tells you the mean, variance, canonical link, and deviance in a single stroke.

The taxonomy organizes the most common exponential-family distributions into continuous (left panel) and discrete (right panel) groups. For each distribution, the canonical parameter θ, cumulant function b(θ), variance function V(μ), and canonical link are displayed. The variance function V(μ) is the key distinguishing feature: it determines how the variance changes with the mean, which is a fundamental assumption one checks against data.

Distribution	θ (canonical)	b(θ)	V(μ)	Canonical Link	φ
Normal	`μ`	`θ²/2`	`1`	Identity: g(μ) = μ	`σ²`
Bernoulli	`log(μ/(1−μ))`	`log(1+eᶿ)`	`μ(1−μ)`	Logit: g(μ) = log(μ/(1−μ))	`1`
Poisson	`log(μ)`	`eᶿ`	`μ`	Log: g(μ) = log(μ)	`1`
Gamma	`−1/μ`	`−log(−θ)`	`μ²`	Reciprocal: g(μ) = 1/μ	`1/α`
Inv. Gaussian	`−1/(2μ²)`	`−√(−2θ)`	`μ³`	Inverse squared: g(μ) = 1/μ²	`λ⁻¹`

Exponential-family reference table for common GLM distributions

An important pattern emerges from this table: the variance function V(μ) ascends in complexity as we move from Normal (V = 1, constant variance) through Poisson (V = μ) and Gamma (V = μ²) to Inverse Gaussian (V = μ³). This sequence corresponds to the Tweedie family V(μ) = μ^p, where p = 0 for Normal, p = 1 for Poisson, p = 2 for Gamma, and p = 3 for Inverse Gaussian. Intermediate values of p (1 < p < 2) yield compound Poisson–Gamma distributions, useful for modeling insurance claims and other data with a point mass at zero plus a continuous positive component.

SECTION 6

Worked Example — Poisson GLM for Count Data

Suppose we are modeling the number of insurance claims Y filed per policyholder per year, with a single covariate x representing the policyholder's age (centered at 40). We assume Y ~ Poisson(μ) and use the canonical log link. We have estimated β₀ = −1.2 and β₁ = 0.03 from a dataset of n = 500 policyholders. Let us trace through the GLM structure and compute the predicted mean and variance for a 55-year-old policyholder.

Poisson GLM: Predicting Insurance Claims

1
Step 1 — Identify the Exponential-Family Components

The Poisson distribution has PMF f(y | μ) = e^−μ μ^y / y!. Rewriting in exponential-family form: f(y | θ) = (1/y!) · exp{y·log(μ) − μ}. Therefore θ = log(μ), b(θ) = e^θ = μ, φ = 1, and a(y, φ) = 1/y!.

θ = log(μ), b(θ) = e^θ, V(μ) = μ
2
Step 2 — Verify the Mean–Variance Relationship

Compute b′(θ) = d(e^θ)/dθ = e^θ = μ. This confirms E[Y] = μ. Then b″(θ) = e^θ = μ, so Var(Y) = φ · b″(θ) = 1 · μ = μ. The Poisson's defining equidispersion property (mean equals variance) falls out automatically.

E[Y] = μ, Var(Y) = μ (equidispersion confirmed)
3
Step 3 — Compute the Linear Predictor for x = 15

For a 55-year-old, the centered covariate is x = 55 − 40 = 15. The linear predictor is η = β₀ + β₁x = −1.2 + 0.03 × 15 = −1.2 + 0.45 = −0.75.

η = −0.75
4
Step 4 — Apply the Inverse Link to Obtain μ̂

The canonical link for Poisson is g(μ) = log(μ), so the inverse link is μ = g⁻¹(η) = exp(η). Thus μ̂ = exp(−0.75) ≈ 0.4724. The model predicts that a 55-year-old policyholder files approximately 0.47 claims per year.

μ̂ = exp(−0.75) ≈ 0.4724 claims/year
5
Step 5 — Compute the Predicted Variance and Interpret β₁

Since Var(Y) = μ for Poisson, the predicted variance is also 0.4724. The coefficient β₁ = 0.03 has a multiplicative interpretation on the mean: a one-year increase in age multiplies the expected claim count by exp(0.03) ≈ 1.030, a 3.0% increase. Equivalently, a 10-year increase in age multiplies the mean by exp(0.30) ≈ 1.350, a 35.0% increase.

Var(Y) ≈ 0.4724; each additional year of age increases expected claims by ≈ 3.0%

SECTION 7

Canonical vs. Non-Canonical Links — Strengths & Limitations

While the canonical link offers compelling theoretical advantages, practical modeling often benefits from alternative link functions. For instance, logistic regression uses the canonical logit link, but the probit link (based on the normal CDF) and the complementary log-log link are common alternatives for binary data, each carrying different assumptions about the latent process. Similarly, for positive continuous data, the gamma distribution is often paired with a log link rather than its canonical reciprocal link because the log link guarantees positive fitted values and has a more natural multiplicative interpretation. The choice of link function reflects both the modeler's beliefs about the relationship between covariates and the response and the practical requirements of the application.

Property	Canonical Link	Non-Canonical Link
Score equations	X⊤(Y − μ̂) = 0; clean, closed-form sufficient statistic X⊤Y	Involve extra weight factor (dμ/dη); no simple sufficient statistic for β
Log-likelihood concavity	Globally concave in β; unique MLE guaranteed	Generally concave in practice but not guaranteed; multiple modes possible in theory
Information matrix	Observed = Expected information; simplifies inference	Observed ≠ Expected; both should be examined for robust inference
Interpretability	Natural for the distribution (e.g., logit = log-odds) but reciprocal link for gamma is less intuitive	May be more interpretable (e.g., log link for gamma gives multiplicative effects)
Boundary constraints	Canonical link always maps to the natural parameter space, so μ respects distributional support	Modeler must verify that g⁻¹(η) stays within the support of μ for all η
Convergence (IRLS)	Typically fast and reliable due to concavity	May require more iterations or step-halving; occasional non-convergence

Comparison of canonical and non-canonical link functions in GLMs

⚖️

PRACTICAL GUIDANCE

Think of the canonical link as the 'default setting' of a sophisticated instrument—engineered to perform optimally under standard conditions. Just as an experienced photographer sometimes departs from the camera's default exposure to handle unusual lighting, an experienced statistician may choose a non-canonical link when the data or scientific context warrant it. The canonical link is always a sensible starting point, but model fit, interpretability, and domain knowledge should guide the final choice. Use residual diagnostics, AIC, and goodness-of-fit tests to compare competing link specifications empirically.

SECTION 8

Connection to Advanced Theory

The exponential-family GLM framework is not an endpoint but a launching pad for more flexible models. Many advanced methods can be understood as extensions that relax one or more of the GLM assumptions while preserving the core exponential-family structure. Understanding the baseline framework makes these extensions both intelligible and natural.

Standard GLM Feature	Extension	What Changes
Fixed linear predictor η = x⊤β	GAMs (Generalized Additive Models)	Replace x⊤β with Σ sⱼ(xⱼ); smooth nonparametric functions via penalized splines
Independence of observations	GLMMs (Generalized Linear Mixed Models)	Add random effects b ~ N(0, D) to η; handle clustering and repeated measures
Single exponential-family distribution	Quasi-likelihood	Specify only mean and variance functions; no full distributional assumption needed
Fully specified parametric model	GEE (Generalized Estimating Equations)	Marginal model with working correlation structure; robust sandwich standard errors
Unpenalized MLE for β	Regularized GLMs (Lasso, Ridge, Elastic Net)	Add penalty λ\|\|β\|\| to negative log-likelihood; variable selection and shrinkage
Frequentist point estimation	Bayesian GLMs	Place priors on β and φ; posterior inference via MCMC or variational methods; exponential-family conjugacy aids computation

Extensions of the standard GLM framework

A particularly elegant connection arises in Bayesian inference: exponential-family likelihoods admit conjugate priors, meaning that the posterior distribution has the same functional form as the prior. For example, a Poisson likelihood paired with a Gamma prior on μ yields a Gamma posterior. This conjugacy—a direct consequence of the exponential-family structure—makes Bayesian updating analytically tractable in simple settings and accelerates MCMC convergence in complex hierarchical models. The exponential family is thus not merely a classification scheme but the algebraic backbone of modern statistical computation.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain why the Uniform(0, 1) distribution does not belong to the exponential family, whereas the Uniform(0, θ) distribution—despite being a seemingly similar parametric family—also fails to qualify. What specific property of the exponential family is violated in each case?

Show answer

The Uniform(0, 1) distribution has no unknown parameters, so the concept of a natural parameter is vacuous—it is a fixed distribution, not a parametric family. The Uniform(0, θ) distribution does have a parameter θ, but its support {y : 0 ≤ y ≤ θ} depends on θ. A requirement for membership in the exponential family (in its standard one-parameter form) is that the support of the density must not depend on the parameter. Formally, the density f(y | θ) = (1/θ) · I(0 ≤ y ≤ θ) includes an indicator function I(0 ≤ y ≤ θ) that cannot be factored into the exp{yθ − b(θ)} form because it couples y and θ in a non-exponential way. This means the sufficient statistic is the order statistic Y₍ₙ₎ = max(Yᵢ), not a sum Σyᵢ, which is the hallmark of exponential families.

Question 2: BASIC CALCULATION

Write the Bernoulli(p) probability mass function in exponential-family form. Identify θ, b(θ), and φ explicitly. Then verify that b′(θ) = μ = p and b″(θ) = Var(Y)/φ = p(1 − p).

Show answer

The Bernoulli PMF is f(y | p) = p^y(1 − p)^1−y for y ∈ {0, 1}. Taking the exponential: f(y | p) = exp{y·log(p/(1−p)) + log(1−p)}. Thus θ = log(p/(1−p)) (the log-odds), b(θ) = log(1 + e^θ), φ = 1, and a(y, φ) = 1. Computing: b′(θ) = e^θ/(1 + e^θ) = p, confirming E[Y] = p. Then b″(θ) = e^θ/(1 + e^θ)² = p(1 − p), confirming Var(Y) = φ · b″(θ) = p(1 − p).

Question 3: INTERMEDIATE

A researcher fits a Poisson GLM with log link to model the number of species y observed in ecological plots, using elevation x₁ (in km) and rainfall x₂ (in 100 mm/year) as covariates. The fitted model is log(μ̂) = 3.1 − 0.8x₁ + 0.4x₂. (a) Compute the predicted number of species for a plot at 2 km elevation with 5 (×100 mm) annual rainfall. (b) Interpret the coefficient −0.8 in terms of a multiplicative effect. (c) If the observed count for such a plot is y = 12, compute the deviance contribution from this observation.

Show answer

(a) η = 3.1 − 0.8(2) + 0.4(5) = 3.1 − 1.6 + 2.0 = 3.5. Therefore μ̂ = exp(3.5) ≈ 33.12 species. (b) The coefficient −0.8 means that each additional kilometer of elevation multiplies the expected species count by exp(−0.8) ≈ 0.449, a 55.1% decrease, holding rainfall constant. (c) The Poisson deviance contribution for a single observation is 2[y·log(y/μ̂) − (y − μ̂)]. Substituting: 2[12·log(12/33.12) − (12 − 33.12)] = 2[12·(−1.015) − (−21.12)] = 2[−12.18 + 21.12] = 2(8.94) = 17.88. This large deviance contribution indicates that y = 12 is substantially lower than the model predicts, warranting investigation of possible overdispersion or missing covariates.

Question 4: APPLIED

A clinical trial models the probability of adverse events using a logistic regression (Bernoulli GLM with logit link). The model includes treatment group (x₁ = 1 for drug, 0 for placebo) and age (x₂, centered at 60). Estimated coefficients are β₀ = −2.5, β₁ = 0.6, β₂ = 0.04. (a) Compute the predicted probability of an adverse event for a 70-year-old patient on the drug. (b) Compute the odds ratio for drug vs. placebo, adjusted for age. (c) The researchers consider using a probit link instead. What would change in the estimation procedure, and would you expect materially different predictions?

Show answer

(a) For a 70-year-old on the drug: x₁ = 1, x₂ = 10. η = −2.5 + 0.6(1) + 0.04(10) = −2.5 + 0.6 + 0.4 = −1.5. The predicted probability is p̂ = 1/(1 + exp(1.5)) = 1/(1 + 4.4817) ≈ 0.1824, or about 18.2%. (b) The odds ratio for drug vs. placebo is exp(β₁) = exp(0.6) ≈ 1.822, meaning the odds of an adverse event are 82.2% higher for the drug group, adjusted for age. (c) Switching to a probit link replaces the logistic CDF with the standard normal CDF: P(Y = 1) = Φ(η). The IRLS algorithm would use different weights—derived from the normal density φ(η) rather than the logistic variance μ(1−μ)—and the coefficients would be on a different scale (roughly βₚᵣₒᵦᵢₜ ≈ βₗₒ_gᵢₜ/1.7). For moderate probabilities (0.1 to 0.9), logit and probit give nearly identical predicted probabilities, so material differences are unlikely unless the data lie heavily in the tails.

Question 5: CRITICAL THINKING

Consider a GLM with a Gamma response distribution. (a) Derive the canonical link by starting from the exponential-family form of the Gamma density. (b) Explain why, despite the canonical link being g(μ) = 1/μ, the log link is overwhelmingly preferred in practice for Gamma GLMs. (c) When a Gamma GLM with log link is fitted, demonstrate that the score equations no longer simplify to X⊤(Y − μ̂) = 0, and identify the additional factor that appears.

Show answer

(a) The Gamma(α, β) density can be parameterized as f(y | μ, α) = [αᵅ / Γ(α)] ·

y^{α−1}

· exp(−αy/μ), where μ = αβ is the mean. Setting φ = 1/α and θ = −1/μ, we get f(y) = a(y, φ)·exp{[yθ − b(θ)]/φ} with b(θ) = −log(−θ). The canonical link is g(μ) = θ = −1/μ, or equivalently g(μ) = 1/μ (the reciprocal). (b) The reciprocal link is problematic because η = 1/μ means μ = 1/η, which can produce negative fitted values μ̂ when η < 0—a violation of the Gamma's positive support. The log link, g(μ) = log(μ), guarantees μ = exp(η) > 0 for all real η and provides a natural multiplicative interpretation of coefficients. (c) With a non-canonical link, the score equations become X⊤D(Y − μ̂) = 0, where D = diag(dμᵢ/dηᵢ · 1/V(μᵢ)) is a diagonal matrix involving both the derivative of the inverse link and the variance function. For the Gamma with log link, dμ/dη = μ (since μ = exp(η)) and V(μ) = μ², so the weight factor is μᵢ/μᵢ² = 1/μᵢ. The score equations become X⊤diag(1/μ̂ᵢ)(Y − μ̂) = 0, which is no longer the simple X⊤(Y − μ̂) = 0 of the canonical case.

SUMMARY

Summary — Exponential Family & Link Functions

Summary

The exponential family provides a unifying algebraic framework for distributions commonly encountered in statistical modeling, including the Normal, Bernoulli, Poisson, and Gamma. Every member can be written in the form f(y | θ, φ) = a(y, φ) · exp{[yθ − b(θ)] / φ}, from which the mean (μ = b′(θ)) and variance (Var = φ · b″(θ)) are obtained by differentiation of the cumulant function b(θ). This mean–variance relationship, captured by the variance function V(μ), is what distinguishes each distribution and drives the weighting in estimation.

A generalized linear model (GLM) pairs an exponential-family distribution with a link function g(μ) = η that maps the mean onto the linear predictor η = x⊤β. The canonical link (where g(μ) = θ) yields particularly elegant score equations X⊤(Y − μ̂) = 0, ensures global concavity of the log-likelihood, and produces X⊤Y as a sufficient statistic. Non-canonical links sacrifice these properties but may offer better interpretability or respect boundary constraints more naturally. Parameters are estimated via iteratively reweighted least squares (IRLS), a Fisher scoring algorithm that unifies estimation across all GLMs. This framework extends naturally to GAMs, GLMMs, GEEs, regularized models, and Bayesian formulations, making it the central organizing principle of modern regression analysis.