From Entropy to Distributions: The Unified Logic of Exponential, Power-Law, Normal, and Boltzmann

The central claim of this article: the exponential, power-law, normal, and Boltzmann distributions are not four separate distributions, but four manifestations of one and the same principle (maximum entropy) under different constraints.

\[\text{Maximize } H = -\!\int p(x)\ln p(x)\,dx \quad \text{subject to constraints} \quad \Longrightarrow \quad p(x) \propto e^{-\lambda\, f(x)}\]

Substitute four different \(f(x)\) into that single equation, and the four distributions drop out one after another:

Distribution	Constrained quantity \(f(x)\)	Analytic form (\(e^{-\lambda f(x)}\) \(\xrightarrow{\text{normalize}}\) \(p(x)\))
Exponential	\(x\) — the mean (\(x \ge 0\))	\(e^{-\lambda x} \;\to\; \lambda e^{-\lambda x}\)
Power law / Pareto	\(\ln x\) — the geometric mean (\(x \ge x_{\min}\))	\(e^{-\lambda \ln x} = x^{-\lambda} \;\to\; (\lambda-1)\,x_{\min}^{\lambda-1}\,x^{-\lambda}\)
Normal	\(x,\ x^2\) — mean and variance	\(e^{-\lambda_1 x-\lambda_2 x^2} \;\to\; \dfrac{1}{\sqrt{2\pi}\,\sigma}\,e^{-(x-\mu)^2/2\sigma^2}\)
Boltzmann	\(E(x)\) — energy	\(e^{-\lambda E(x)} \;\to\; \dfrac{1}{Z}\,e^{-E(x)/k_BT}\)

The second row is worth a stare: \(e^{-\lambda\ln x}\) is \(x^{-\lambda}\)—switch the constraint from “the mean of \(x\)” to “the mean of \(\ln x\)”, and the exponential turns into a power law on the spot. Likewise, adding an \(x^2\) term in the third row completes the square into a Gaussian, and reading \(\lambda\) as \(1/k_BT\) in the fourth row gives the physicist’s Boltzmann factor. There is no leap between the four distributions—only one substitution of \(f(x)\).

---

1 Why maximum entropy?

Before maximizing entropy, answer a more fundamental question first: why should we maximize entropy at all? Two independent roads lead to it—the physical counting argument, and the statistical unbiased-inference argument.

1.1 The counting argument: entropy = log of the number of arrangements

Make \(N\) observations, each landing in one of \(k\) categories. Let \(n_i\) be the number of times category \(i\) occurs (\(i = 1, \dots, k\), with \(n_1 + n_2 + \cdots + n_k = N\)). The macrostate is the count vector \((n_1, \dots, n_k)\)—recording only “how many times each category occurs”; the microstate is the specific observation sequence—recording additionally “which observation was which category.” The number of microscopic arrangements corresponding to one macrostate is the multinomial coefficient:

\[W = \frac{N!}{n_1!\,n_2!\cdots n_k!}\]

Coin tossing is the special case \(k=2\): only heads/tails, \(n_1 = k_{\text{heads}}\), \(n_2 = N - k_{\text{heads}}\), so

\[W = \frac{N!}{k_{\text{heads}}!\,(N-k_{\text{heads}})!} = \binom{N}{k_{\text{heads}}}\]

Proof that \(\ln W \approx N \cdot H\):

\[ \begin{aligned} \ln W &= \ln N! - \sum_i \ln n_i! \\[2pt] &\approx \bigl(N\ln N - N\bigr) - \sum_i \bigl(n_i \ln n_i - n_i\bigr) && \text{Stirling: } \ln n! \approx n\ln n - n \\[2pt] &= N\ln N - \sum_i n_i \ln n_i && \textstyle\sum_i n_i = N \text{, the two linear terms cancel} \\[2pt] &= \sum_i n_i \ln N - \sum_i n_i \ln n_i && \textstyle N\ln N = \bigl(\sum_i n_i\bigr)\ln N \text{, constant } \ln N \text{ pulled into the sum} \\[2pt] &= \sum_i n_i \bigl(\ln N - \ln n_i\bigr) && \text{same summation index, combine and factor out } n_i \\[2pt] &= -\sum_i n_i \ln\frac{n_i}{N} && \textstyle\ln N - \ln n_i = -\ln\frac{n_i}{N} \text{ (quotient rule for logs), minus sign pulled out} \\[2pt] &= N \cdot \left(-\sum_i p_i \ln p_i\right) && p_i = n_i / N \\[2pt] &= N \cdot H(p) \qquad \blacksquare \end{aligned} \]

That is, \(W \approx e^{N H}\): entropy is the “arrangement-count exponent” amortized over each observation. Even if two distributions differ in entropy by only a tiny amount, their arrangement counts differ by a factor of \(e^{N\cdot\Delta H}\)—an astronomical number when \(N\) is large. So “the maximum-entropy distribution” is not an abstract preference but the overwhelmingly more probable macrostate.

## ── Counting intuition: higher entropy = more microscopic arrangements = overwhelmingly more likely ──

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))

# Left: 30 coin flips; arrangements W = C(30, k) for k heads
n_flip <- 30
k <- 0:n_flip
W <- choose(n_flip, k)
plot(k / n_flip, W, type = "h", lwd = 3, col = "steelblue",
     main = paste0("W = C(", n_flip, ", k): microstates per macrostate"),
     xlab = "Fraction of heads  k/N", ylab = "Number of arrangements W",
     cex.lab = 1.35, cex.main = 1.15)
points(0.5, choose(n_flip, n_flip/2), pch = 16, col = "tomato", cex = 1.4)
text(0.02, 1.3e8, adj = 0,
     paste0("k/N = 1/2:\nW = ", format(choose(n_flip, n_flip/2), big.mark = ",")),
     col = "tomato", cex = 1.15)
text(1.02, 3.2e7, "k/N = 1:\nW = 1\n(bar invisible)", col = "gray40", cex = 1.05, adj = 1)
arrows(0.97, 2.0e7, 1.0, 2e6, length = 0.08, col = "gray40", lwd = 1.5)

# Right: (1/N) ln C(N, pN) converges to the binary entropy H(p) as N grows
p_grid <- seq(0.001, 0.999, length.out = 400)
Ns   <- c(10, 50, 1000)
cols <- c("forestgreen", "purple", "steelblue")
plot(NA, xlim = c(0, 1), ylim = c(0, 0.75),
     main = "ln(W)/N converges to entropy H(p)",
     xlab = "p (fraction of heads)", ylab = "ln(W) / N",
     cex.lab = 1.35, cex.main = 1.15)
for (i in seq_along(Ns)) {
  N_i <- Ns[i]
  lines(p_grid, lchoose(N_i, round(p_grid * N_i)) / N_i,
        col = cols[i], lwd = 2)
}
H_bin <- -p_grid * log(p_grid) - (1 - p_grid) * log(1 - p_grid)
lines(p_grid, H_bin, col = "tomato", lwd = 3, lty = 2)
legend("bottom", c(paste0("N = ", Ns), "H(p) = -p ln p - (1-p) ln(1-p)"),
       col = c(cols, "tomato"), lwd = c(2, 2, 2, 3),
       lty = c(1, 1, 1, 2), bty = "n", cex = 0.95)

Left: each bar is a macrostate (fraction of heads \(k/N\)), and its height is the corresponding number of microscopic arrangements \(W = \binom{30}{k}\). Note that the bars at \(k/N = 0\) and \(1\) are drawn but invisible—their height is \(W = 1\) (all-tails/all-heads has only one arrangement), which is visually zero against a vertical axis reaching 155 million. That is exactly the message: the arrangement count at \(k/N = 1/2\) is 150 million times that at the extremes; when you “see” half heads from a fair coin, it is not that some force pushes it there—it is pure arrangement-count dominance. Right: \(\frac{1}{N}\ln W\) converges to \(H(p)\) as \(N\) grows, the numerical verification of the Stirling proof above.

1.2 How wide is the peak? Second-order expansion and the \(1/\sqrt{N}\) window

There is something in the figure above that invites suspicion: for all the talk of “dominance,” why do the middle bars differ so little? For instance \(\binom{30}{14}/\binom{30}{15} = 0.94\), only a 6% difference. This is no contradiction—it is a necessary consequence of the second-order expansion.

\(W_1/W_2 = e^{N\Delta H}\) really is the ratio of the \(W\)’s themselves, but what matters is how large \(\Delta H\) is. \(H(p)\) has a smooth maximum at \(p = 1/2\), where the first derivative vanishes, so the Taylor expansion starts at second order (write \(\delta = p - 1/2\)):

\[\Delta H = H(1/2) - H(1/2+\delta) \approx \tfrac{1}{2}\bigl|H''(1/2)\bigr|\,\delta^2 = 2\delta^2\]

So near the peak:

\[\frac{W(p)}{W_{\max}} \approx e^{-2N\delta^2}\]

Two conclusions read off directly:

• Adjacent bars differ little: at \(\delta = 1/30\), \(e^{-2N\delta^2} \approx e^{-1/15} \approx 0.94\)—the “flat middle” you see is the signature of the vanishing first derivative;
• Window width \(\sim 1/\sqrt{N}\): for \(e^{-2N\delta^2}\) to fall to negligible values requires \(\delta \gtrsim 1/\sqrt{N}\). This is exactly the CLT’s \(\sqrt{N}\) fluctuation—arrangement counting and the central limit theorem are the same thing here.

So “dominance” is about a fixed macroscopic difference (fixed \(\Delta H\)) as \(N\) grows; at fixed \(N\) there is always an “undominated” window of width \(O(1/\sqrt{N})\) near the peak, and as \(N\to\infty\) the window itself shrinks to zero.

## ── Dominance as a function of N: window shrinks + fixed ΔH collapses exponentially ──

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))

# Left: W/W_max sharpens into a spike as N grows (width ~ 1/sqrt(N))
Ns_w  <- c(30, 300, 3000)
cols_w <- c("steelblue", "purple", "tomato")
plot(NA, xlim = c(0.2, 0.8), ylim = c(0, 1.08),
     main = "Relative arrangements W / W_max",
     xlab = "p = k/N", ylab = "W(k) / W(N/2)",
     cex.lab = 1.35, cex.main = 1.15)
for (i in seq_along(Ns_w)) {
  N_i <- Ns_w[i]
  k_i <- 0:N_i
  rel <- exp(lchoose(N_i, k_i) - lchoose(N_i, round(N_i / 2)))
  lines(k_i / N_i, rel, col = cols_w[i], lwd = 2.5)
}
abline(v = 0.5, col = "gray70", lty = 3)
legend("topleft", paste0("N = ", Ns_w),
       col = cols_w, lwd = 2.5, bty = "n", cex = 1.1,
       title = "width ~ 1/sqrt(N)", title.col = "gray30")

# Right: fixed macroscopic contrast p=0.6 vs p=0.5, W ratio collapses exponentially in N
Ns_r <- seq(10, 1000, by = 10)
ratio <- sapply(Ns_r, function(N_i)
  exp(lchoose(N_i, round(0.6 * N_i)) - lchoose(N_i, round(0.5 * N_i))))
plot(Ns_r, ratio, type = "l", lwd = 2.5, col = "steelblue", log = "y",
     main = "Fixed contrast: W(p=0.6) / W(p=0.5)",
     xlab = "N", ylab = "Ratio (log scale)",
     cex.lab = 1.35, cex.main = 1.15)
dH <- (-0.5*log(0.5)-0.5*log(0.5)) - (-0.6*log(0.6)-0.4*log(0.4))
lines(Ns_r, exp(-Ns_r * dH), col = "tomato", lwd = 2, lty = 2)
pts_N <- c(30, 1000)
pts_r <- sapply(pts_N, function(N_i)
  exp(lchoose(N_i, round(0.6 * N_i)) - lchoose(N_i, round(0.5 * N_i))))
points(pts_N, pts_r, pch = 16, col = "forestgreen", cex = 1.5)
text(pts_N[1], pts_r[1], pos = 4,
     paste0("N = 30:  ratio = ", round(pts_r[1], 2)), col = "forestgreen", cex = 1.15)
text(pts_N[2], pts_r[2] * 1e3, pos = 2,
     paste0("N = 1000:  ratio = ", format(pts_r[2], digits = 2)),
     col = "forestgreen", cex = 1.15)
legend("bottomleft", c("exact (binomial)", "e^(-N dH),  dH = 0.0201"),
       col = c("steelblue", "tomato"), lwd = c(2.5, 2), lty = c(1, 2),
       bty = "n", cex = 1.0)

Left: dividing each \(N\)’s arrangement count by its own peak, \(N = 30 \to 3000\), the window shrinks as \(1/\sqrt{N}\). Right: the fixed comparison \(p = 0.6\) vs \(p = 0.5\) (\(\Delta H \approx 0.0201\)) gives a straight line \(e^{-N\Delta H}\) on the log axis—at \(N=30\) the difference is only a factor of two (0.56, unremarkable), but at \(N=1000\) it is already \(10^{-9}\). The exact binomial computation (blue) and \(e^{-N\Delta H}\) (red dashed) nearly coincide.

1.3 Historical note: the same \(H\), discovered twice

The derivation above looks like information theory, but historically physics got there first—the quantity \(-\sum p_i \ln p_i\) was discovered independently twice, by completely different routes:

• 1730–1733, de Moivre and Stirling: de Moivre, studying the binomial distribution of coin tosses, first derived the shape \(n! \sim C\, n^{n+1/2} e^{-n}\); Stirling pinned down the constant \(C = \sqrt{2\pi}\) in Methodus Differentialis (1730), and de Moivre then applied the tool to the normal approximation of the binomial in a 1733 pamphlet. The tool was in place. (Moivre 1733; Stirling 1730)

• 1877, Boltzmann: the proof in this section is the work of his paper that year (On the Relationship between the Second Law of Thermodynamics and Probability Theory)—chop the energy of gas molecules into discrete quanta, count the arrangements of each allocation (his Komplexionen, i.e. \(W\)), apply Stirling, maximize under constraints → the most probable macrostate. \(S = k_B \ln W\) was later carved on his tombstone; divide by \(N\) and you get \(-\sum p_i \ln p_i\). (Boltzmann 1877)

  The mathematics looks like an undergraduate exercise today, but the difficulty of a formula ≠ the difficulty of the insight—every step was heresy at the time: (a) claiming “entropy is a probabilistic thing and the second law is merely a statistical regularity” amounts to saying entropy increase can be violated (just with vanishing probability), which drew fierce attacks (Loschmidt’s reversibility paradox, Zermelo’s recurrence paradox—a twenty-year controversy); (b) the whole argument rests on the reality of atoms, which the mainstream led by Mach and Ostwald did not accept; Boltzmann fought this battle for the second half of his life and took his own life in 1906, not living to see Perrin’s Brownian-motion experiments confirm atomism a few years later; (c) the energy discretization \(\epsilon, 2\epsilon, 3\epsilon,\dots\) was to him just a computational trick to be taken to a limit at the end, but Planck, treating black-body radiation in 1900, found that not taking the limit was the right move—quantum theory grew out of this “trick.” The notation \(S = k\ln W\) and the constant \(k\) were in fact both given by Planck, who named \(k\) the Boltzmann constant. (Cercignani 2006; Planck 1901)

• 1948, Shannon (A Mathematical Theory of Communication): took the axiomatic route—requiring the uncertainty measure to be continuous in \(p\), to increase monotonically with the number of categories when they are equiprobable, and to decompose additively under step-by-step choices, he proved the unique solution is \(-K\sum p_i \log p_i\). No arrangement counts, no Stirling in the derivation. But the typical sequences in his paper loop back to counting: among sequences of length \(N\), about \(e^{NH}\) are “typical”—exactly the information-theoretic incarnation of \(W \approx e^{NH}\). (Shannon 1948)

• Naming anecdote: a commonly repeated version says that von Neumann suggested Shannon call the quantity “entropy,” because “no one really understands what entropy is.” The story is traceable to a 1971 retelling rather than a verifiable conversation record, so it is best treated as a second-hand anecdote. (Tribus and McIrvine 1971)

• 1957, Jaynes (Information Theory and Statistical Mechanics): formally welded the two roads together—the maximum entropy of statistical mechanics is not a physical law but a special case of an inference principle. That is the next subsection. (Jaynes 1957)

1.4 Jaynes: maximum entropy = the most honest inference

The second road speaks not of physics but of inference. (Jaynes 1957) When you know only partial information about a distribution (say, only the mean) but must choose a complete distribution:

• Choosing the one with maximum entropy = using only the known constraints, smuggling in no extra assumptions;
• Choosing any other distribution = secretly claiming to know something you actually do not (e.g. “the distribution is bimodal”).

So maximum entropy is not nature’s mysterious preference but the unique consistent way to stay maximally ignorant beyond the constraints. Physical systems “choose” it because of arrangement-count dominance; statisticians “choose” it out of honesty—the two roads converge, both pointing to the same variational problem:

\[\max_p \; H[p] \quad \text{s.t. known constraints} \;\Longrightarrow\; p(x) \propto e^{-\lambda f(x)}\]

Chapters 3 through 6 that follow are just plugging in different \(f(x)\).

Watch next (optional). Reinventing Entropy | Compression is Intelligence Part 1 (3Blue1Brown) reconstructs intuition for entropy from compression; it is broader than this article’s maximum-entropy thread and works best after this section.

---

2 Entropy: why these distributions and not others?

2.1 Definition of information entropy

For a continuous distribution \(p(x)\), the differential entropy is defined as:

\[H[p] = -\int p(x) \ln p(x)\, dx\]

The larger \(H\), the more “spread out,” more “uncertain,” and more “unbiased” the distribution.

(There is in fact a “ruler” hidden in this definition: spread out relative to what? In the continuous case, changing units changes the numerical value of \(H\). That issue only needs facing head-on in the closing chapter, First, the ruler; until then the conventional reading is fine.)

A note on notation: why information entropy is written \(H\) and not \(S\). Shannon borrowed the letter from Boltzmann’s H-theorem—he says so explicitly in the 1948 paper: this \(H\) is the \(H\) in \(H = \int f\ln f\,dv\), the quantity Boltzmann proved to be monotonically decreasing (so \(S = -k_B H\); Shannon kept the letter and flipped the sign). As for why Boltzmann used H, one account holds that it was originally a capital Greek Eta (Η)—entropy is written \(\eta\) in some traditions, and a capital Eta is typographically indistinguishable from a Latin H.

What matters is that the two are not the same quantity; they differ by \(k_B\) and by dimension:

	Symbol	Definition	Dimension
Information entropy (used throughout)	\(H\)	\(-\sum p\ln p\)	dimensionless (nats; bits if \(\log_2\))
Thermodynamic entropy	\(S\)	\(-k_B\sum p\ln p\)	J/K

\[S = k_B H\]

This document follows the conventions of each field: \(H\) when discussing inference, \(S\) when discussing thermodynamics (temperature, the Boltzmann distribution), converting between them with the identity above. Physicists routinely write \(S\) for this very same functional as well (the Gibbs entropy, Tsallis entropy \(S_q\)), so \(H\) versus \(S\) is a difference of convention, not of correctness. One more caution: the \(H\) of thermochemistry is enthalpy, an unrelated quantity that never appears here.

2.2 Reading the integrand: \(p(x) \times (-\ln p(x))\)

The product \(-p\ln p\) in the definition is worth taking apart. Rewrite it as \(H = \int p(x)\cdot\bigl(-\ln p(x)\bigr)\,dx = E\bigl[-\ln p(X)\bigr]\):

• \(-\ln p(x)\) is the information content (surprise) at \(x\): the smaller \(p\), the more surprising, diverging to \(+\infty\) in the tails;
• \(p(x)\) is the probability weight: the frequency with which that surprise actually occurs;
• entropy = average surprise.

## ── Decompose the integrand: p(x) × (-ln p(x)) ──────────────────────────

x_e  <- seq(-4, 4, length.out = 2000)
p_e  <- dnorm(x_e)          # probability weight
s_e  <- -log(p_e)           # information content (surprise)
g_e  <- p_e * s_e           # entropy density = product of the two

pts      <- c(0, 1.5, 3)    # three tracked points: center / mid / tail
pt_cols  <- c("tomato", "purple", "forestgreen")

par(mfrow = c(1, 3), mar = c(4.6, 4.8, 3.5, 1))

# Panel 1: probability weight p(x)
plot(x_e, p_e, type = "l", lwd = 3, col = "steelblue",
     main = "Weight  p(x)", xlab = "x", ylab = "p(x)",
     cex.lab = 1.4, cex.main = 1.3)
points(pts, dnorm(pts), pch = 16, col = pt_cols, cex = 1.6)
text(3, dnorm(3) + 0.05, "p almost 0", col = "forestgreen", pos = 2, cex = 1.2)

# Panel 2: information content -ln p(x)
plot(x_e, s_e, type = "l", lwd = 3, col = "steelblue",
     main = "Surprise  -ln p(x)", xlab = "x", ylab = "-ln p(x)",
     cex.lab = 1.4, cex.main = 1.3)
points(pts, -log(dnorm(pts)), pch = 16, col = pt_cols, cex = 1.6)
text(3, -log(dnorm(3)) - 0.6, "huge surprise", col = "forestgreen",
     pos = 2, cex = 1.2)

# Panel 3: product = entropy density, shaded area = H
plot(x_e, g_e, type = "l", lwd = 3, col = "steelblue",
     main = "Product  p(x) * (-ln p(x))", xlab = "x",
     ylab = "p(x) * (-ln p(x))", cex.lab = 1.4, cex.main = 1.3)
polygon(c(x_e, rev(x_e)), c(g_e, rep(0, length(g_e))),
        col = adjustcolor("steelblue", 0.25), border = NA)
points(pts, dnorm(pts) * (-log(dnorm(pts))), pch = 16, col = pt_cols, cex = 1.6)
H_val <- 0.5 * log(2 * pi * exp(1))
text(0, 0.18, paste0("Area = H = ", round(H_val, 3)),
    col = "steelblue", cex = 1.3)
text(3, dnorm(3) * (-log(dnorm(3))) + 0.04,
     "tiny x huge = ~0", col = "forestgreen", pos = 2, cex = 1.2)

Tracking three points (standard normal):

• red point \(x=0\): maximum weight (0.399) but minimum surprise (0.92)—a moderately large product;
• purple point \(x=1.5\): both weight and surprise are moderate—the product is still substantial;
• green point \(x=3\): large surprise (5.4) but nearly zero weight (0.004)—product \(\approx 0\).

The tail does not blow up because \(\lim_{p\to 0} p\ln p = 0\): \(p\) goes to 0 (exponentially fast) always faster than \(\ln p\) diverges (logarithmically). So the contribution to entropy comes mainly from the “moderate-probability” region, and the shaded area in the third panel is exactly \(H = \tfrac{1}{2}\ln(2\pi e) \approx 1.419\).

2.3 Entropy grows as the distribution spreads

## ── Intuition: entropy measures how "spread out" a distribution is ──────────

# Helper: numerically compute differential entropy
h_diff <- function(p_vals, dx) {
  p_vals <- pmax(p_vals, 1e-300)
  -sum(p_vals * log(p_vals)) * dx
}

x_grid <- seq(-6, 6, length.out = 5000)
dx     <- diff(x_grid)[1]

# Normal distributions with different sigma
sigmas  <- c(0.5, 1, 2, 3)
h_vals  <- sapply(sigmas, function(s) h_diff(dnorm(x_grid, 0, s), dx))
h_theory <- 0.5 * log(2 * pi * exp(1) * sigmas^2)

par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))

# Left: wider distribution -> higher entropy
plot(NA, xlim = c(-6, 6), ylim = c(0, 0.85),
     main = "Wider distribution = higher entropy",
     xlab = "x", ylab = "Density")
cols <- c("steelblue", "tomato", "forestgreen", "purple")
for (i in seq_along(sigmas)) {
  lines(x_grid, dnorm(x_grid, 0, sigmas[i]), col = cols[i], lwd = 2.5)
}
legend("topright",
       as.expression(lapply(seq_along(sigmas), function(i)
         bquote(list(sigma == .(sigmas[i]), H == .(round(h_theory[i], 2)))))),
       col = cols, lwd = 2.5, bty = "n", cex = 0.95)

# Right: sigma vs entropy
plot(sigmas, h_theory, type = "b", pch = 16, col = "steelblue", lwd = 2,
     cex = 1.5, main = "Entropy grows with spread",
     xlab = expression(sigma), ylab = "Differential entropy H")
lines(seq(0.3, 3.2, 0.01),
      0.5 * log(2*pi*exp(1)*seq(0.3,3.2,0.01)^2),
      col = "tomato", lwd = 2, lty = 2)
legend("bottomright",
       expression("Numerical",
                  "Analytic:  " * H == frac(1, 2) ~ ln(2 * pi * e * sigma^2)),
       col = c("steelblue","tomato"), lwd = 2, pch = c(16, NA),
       lty = c(1, 2), bty = "n", y.intersp = 1.4)

Intuition summary:

• concentrate all probability at one point → \(H \to -\infty\) (completely certain, lowest entropy)
• spread out uniformly → \(H\) is maximal (completely uncertain, highest entropy)
• under constraints the distribution cannot be perfectly uniform, but it can spread out as much as the constraints allow

2.4 The maximum entropy principle

When you know only certain statistics about the data (mean, variance, …), the most honest thing to do is to choose the distribution that satisfies these constraints but “assumes” as little as possible about everything else—i.e. the distribution with maximum entropy.

(In English, Maximum Entropy, abbreviated MaxEnt, which is what the English labels inside the figures use.)

This is not an arbitrary aesthetic preference but a logical necessity:

• if you choose a lower-entropy distribution, you are secretly assuming information you do not actually have
• the maximum-entropy distribution = the distribution that injects no extra information beyond the constraints

2.5 Lagrange multipliers: from scratch

This section and the next (the geometry and derivation of Lagrange multipliers) are pure methodological groundwork. If you have already learned Lagrange multipliers, or don’t want to see the derivation for now, feel free to skip both sections and go straight to “Constraints determine the distribution” for the conclusion, or even jump directly to the concrete distributions that follow (starting with the “Exponential distribution” chapter)—it won’t affect the later reading.

The next subsection uses Lagrange multipliers, so let’s first make the tool itself clear. The core is a single figure, and the math is just one line.

The shape of the problem. An ordinary extremum problem: where is \(f(x,y)\) largest?—derivative zero, \(\nabla f = 0\). A constrained extremum problem: subject to \(g(x,y) = c\), where is \(f\) largest? You can only walk along the constraint curve, and at the optimum the derivative of \(f\) is generally not zero—you are simply blocked by the constraint.

Geometric intuition (the soul of the whole method). Walk along the constraint curve, watching \(f\):

• if \(\nabla f\) (the direction of steepest increase of \(f\)) has a component along the curve’s tangent, then stepping along the curve improves \(f\)—not yet optimal;
• at the optimum, \(\nabla f\) must be perpendicular to the constraint curve (zero tangential component, nowhere to go).

And \(\nabla g\) is always perpendicular to the curve \(g = c\) (along it \(g\) does not change). Two vectors perpendicular to the same curve → parallel:

\[\boxed{\nabla f = \lambda\, \nabla g}\]

\(\lambda\) is the Lagrange multiplier—its entire meaning is the proportionality factor of “the two gradients are parallel at the optimum.”

The operating recipe. Construct the Lagrangian to absorb the constraint into the objective:

\[L(x, y, \lambda) = f(x,y) - \lambda\,\bigl(g(x,y) - c\bigr)\]

Take partial derivatives with respect to all variables (including \(\lambda\)) and set them to zero: the first two equations give \(\nabla f = \lambda\nabla g\), and \(\partial L/\partial\lambda = 0\) automatically restores the constraint. The constrained problem becomes an unconstrained one—that is why it is so useful.

Example: rectangle with perimeter 20, maximum area? Maximize \(f = xy\) subject to \(x + y = 10\):

\[L = xy - \lambda(x + y - 10), \qquad \frac{\partial L}{\partial x} = y - \lambda = 0, \quad \frac{\partial L}{\partial y} = x - \lambda = 0 \;\Rightarrow\; x = y = \lambda = 5\]

A square, area 25.

The figure below draws this as a mountain: the surface is the full \(f = xy\), the whole mountain rising toward the \((10,10)\) direction. The blue rings are contour lines (topographic-map style, \(f\) constant along each); the direction of \(\nabla f\) can be read off at a glance—crossing the contours perpendicularly, pointing toward the next-higher ring, which is the direction of steepest ascent. The constraint line, lifted onto the surface, is a path across the hillside (red); constrained optimization = walking only along this path to find its highest point. Look at the key purple point \((5,5)\): \(\nabla f\) is not zero here—the green arrow, perpendicular to the contours, points uphill, and \(f\) could still grow, but that direction leaves the constraint line and is forbidden; what is zero is only the tangential-along-the-path component of \(\nabla f\) (right panel: the slope of the path’s height profile is zero here). And at the purple point the constraint line is exactly tangent to the contour \(f = 25\)—at the instant the path runs along a contour, \(f\) neither increases nor decreases, precisely the extremum.

Two more things are drawn. The translucent vertical wall is the image of the constraint \(g = 10\) in 3D: \(g(x,y) = x+y\) contains no \(z\), so \(\{x + y = 10\}\) is in 3D a wall extending freely in \(z\) (drawn here only up to the mountain surface), and the red path is the cut where the wall slices into the mountain; the dashed line on the ground is the wall’s base, i.e. the constraint set itself. The orange arrow is \(\nabla g\) lifted into 3D as \((1,1,0)\): it is the wall’s unique normal—there are infinitely many directions perpendicular to a curve in 3D (an entire normal plane), but only one perpendicular to a surface, so “\(\nabla g\) is perpendicular to the constraint” is least likely to be misread when thought of via the wall. The 2D version on the ground says the same thing: \(\nabla g = (1,1)\) has only two components, cannot leave the ground, is perpendicular to the dashed line in the plane, and its direction is likewise unique.

How many dimensions is the gradient? The number of components = the number of inputs. \(f(x,y) = xy\) has two inputs, so \(\nabla f = (y,\ x)\) is a 2D vector, \(\nabla f(5,5) = (5,\ 5)\); \(z = 25\) is the output and takes no component. The graph point \((5,5,25)\) lives in 3D, but \(f\) and its gradient both live in the 2D domain (the “map”). The direction \((1,1)\) gives steepest ascent, and the magnitude gives the slope: \(|\nabla f| = \sqrt{50} \approx 7.07\), i.e. walking 1 horizontally along the steepest direction, \(f\) rises by 7.07.

Let’s verify by hand. Walking in the \((1,1)\) direction, first normalize to the unit vector \(\bigl(\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2}\bigr)\) (the vector \((1,1)\) has length \(\sqrt2\); without normalizing, \(t\) would not be the distance traveled); after traveling distance \(t\) the position is \(x(t) = y(t) = 5 + \tfrac{t}{\sqrt2}\), so:

\[f(t) = x(t)\,y(t) = \Bigl(5 + \tfrac{t}{\sqrt2}\Bigr)^{2} = 25 + \tfrac{10}{\sqrt2}\,t + \tfrac{t^2}{2} \;\Rightarrow\; f'(0) = \tfrac{10}{\sqrt2} = \sqrt{50} \approx 7.07\ \checkmark\]

The same result follows in one step from the directional-derivative formula: \(\nabla f \cdot \hat u = (5,5)\cdot\bigl(\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2}\bigr) = \tfrac{10}{\sqrt2}\)—the dot product is largest when \(\hat u\) is aligned with \(\nabla f\), equal exactly to \(|\nabla f|\). “Gradient magnitude = steepest slope” is not a definition but a computed fact.

Calibrate the legend accordingly: the green arrow is a visualization vector hugging the hillside, and only its horizontal projection is \(\nabla f\); likewise \(\nabla g\) is drawn as \((1,1,0)\) (padding a 0 into the figure). The genuinely 3D gradient belongs to a different function—for \(F(x,y,z) = z - xy\), \(\nabla F = (-y, -x, 1)\), which is perpendicular to the hillside surface itself. General rule: the gradient of an \(n\)-variable function is an \(n\)-dimensional vector, perpendicular to its \((n-1)\)-dimensional level set. The Lagrange condition uses \(\nabla f\).

## ── Geometry of Lagrange multipliers: the constraint line is a "mountain path" on the surface ──

par(mfrow = c(1, 2), mar = c(1, 1, 3.5, 0.5))

# Left: surface z = xy; constraint x + y = 10 lifted onto the surface is a mountain path
x_s <- seq(0, 10, length.out = 41)
y_s <- seq(0, 10, length.out = 41)
z_s <- outer(x_s, y_s)

pm <- persp(x_s, y_s, z_s, zlim = c(0, 100),
            theta = -63, phi = 24, expand = 0.65,
            col = "gray92", border = "gray78",
            xlab = "x", ylab = "y", zlab = "f = xy",
            main = "Full surface f = xy;\nconstraint = a path across the hillside",
            cex.main = 1.1)

# Complete the bounding cube: persp only draws the hidden edges near the origin; fill in the x=10 / y=10 sides
bx <- function(x0, y0, z0, x1, y1, z1)
  lines(trans3d(c(x0, x1), c(y0, y1), c(z0, z1), pm),
        col = "gray20", lty = 3, lwd = 1.8)
bx(10, 10, 0,   10, 10, 100)   # vertical edge at the (10,10) corner
bx(10,  0, 0,   10, 10, 0)     # back edge of the base, x = 10
bx( 0, 10, 0,   10, 10, 0)     # back edge of the base, y = 10
bx(10,  0, 100, 10, 10, 100)   # back edge of the top, x = 10
bx( 0, 10, 100, 10, 10, 100)   # back edge of the top, y = 10

# Vertical wall: g = 10 extends along z (x+y=10, z free); the wall's cut into the mountain = red path
xw <- seq(0, 10, length.out = 120)
wall <- trans3d(c(xw, rev(xw)), c(10 - xw, rev(10 - xw)),
                c(rep(0, 120), rev(xw * (10 - xw))), pm)
polygon(wall$x, wall$y, col = adjustcolor("slateblue", 0.28), border = NA)
# grad g = (1,1,0): the wall's unique normal (horizontal, perpendicular to the wall)
gw <- trans3d(c(7.5, 8.7), c(2.5, 3.7), c(0, 0), pm)
arrows(gw$x[1], gw$y[1], gw$x[2], gw$y[2],
       col = "darkorange", lwd = 3, length = 0.12)
points(trans3d(7.5, 2.5, 0, pm), pch = 16, col = "darkorange", cex = 1.1)
text(trans3d(11.4, 3.0, 12, pm), "grad g = (1,1,0):\nthe wall's\nunique normal",
     col = "darkorange3", cex = 1.0)

# Contour lines on the mountain surface (topographic rings): grad f is everywhere perpendicular to them
cls <- contourLines(x_s, y_s, z_s, levels = c(9, 16, 25, 36, 49, 64, 81))
for (cl in cls) lines(trans3d(cl$x, cl$y, cl$level, pm), col = "steelblue", lwd = 1.4)

x_p <- seq(0, 10, length.out = 200)
# Constraint line (on the ground)
lines(trans3d(x_p, 10 - x_p, 0, pm), col = "gray30", lwd = 2, lty = 2)
# Mountain path: constraint line lifted onto the surface, height f = x(10-x)
lines(trans3d(x_p, 10 - x_p, x_p * (10 - x_p), pm), col = "tomato", lwd = 3.5)
# Non-optimal point (2,8) and optimal point (5,5)
points(trans3d(2, 8, 16, pm), pch = 16, col = "tomato", cex = 1.3)
points(trans3d(5, 5, 25, pm), pch = 16, col = "purple", cex = 1.7)
# Vertical drop line at the optimum
lines(trans3d(c(5, 5), c(5, 5), c(0, 25), pm), col = "purple", lty = 3, lwd = 1.5)
# grad f at (5,5) is nonzero: points uphill (leaving the constraint line)
gr <- trans3d(c(5, 6.6), c(5, 6.6), c(25, 6.6^2), pm)
arrows(gr$x[1], gr$y[1], gr$x[2], gr$y[2], col = "forestgreen", lwd = 3, length = 0.12)
text(trans3d(5.6, 5.6, 76, pm),
     "grad f != 0:\nsteepest ascent,\nperpendicular to contours,\nbut leaves the constraint",
     col = "forestgreen", cex = 1.0)
text(trans3d(2.2, 2.2, 16, pm), "top of the path\n(5, 5),  f = 25", col = "purple", cex = 1.1)
text(trans3d(1.0, 9.4, 30, pm), "still climbing", col = "tomato", cex = 1.05)
text(trans3d(2.2, 7.8, -8, pm), "x + y = 10", col = "gray30", cex = 1.05)

# Right: walking along the constraint line, f = x(10-x) (side view of the mountain path)
par(mar = c(4.6, 5.0, 3.5, 1))
x_c <- seq(0, 10, length.out = 300)
plot(x_c, x_c * (10 - x_c), type = "l", lwd = 3, col = "tomato",
     main = "Height along the path:  f = x (10 - x)",
     xlab = "x  (position on the line)", ylab = "f = xy",
     cex.lab = 1.35, cex.main = 1.15)
points(5, 25, pch = 16, col = "purple", cex = 1.6)
points(2, 16, pch = 16, col = "tomato", cex = 1.4)
arrows(2.3, 16.8, 3.4, 21, col = "tomato", lwd = 2, length = 0.1)
text(3.0, 12.6, "still climbing", col = "tomato", cex = 1.05)
text(5, 22.2, "top: tangential slope = 0", col = "purple", cex = 1.1)

The meaning of \(\lambda\): the “price” of the constraint. A beautiful theorem: the sensitivity of the optimal value to the constraint is the multiplier itself,

\[\frac{df^*}{dc} = \lambda\]

In the example \(\lambda = 5\): relaxing the half-perimeter from 10 to 11, the maximum area rises by about 5 (exact value \(5.5^2 - 5^2 = 5.25\)). (The numerical value of \(\lambda\) is tied to how the constraint is written: multiply \(g\) by 2 and \(\lambda\) halves, \(df^*/dc\) halving in step, so the theorem is self-consistent under any convention; its dimension is \([f]/[g]\).)

This “price of the constraint” has its own name in every field:

Field	Name of \(\lambda\)	Example
Optimization/math	dual variable, KKT multiplier	solution of the LP dual
Analytical mechanics	constraint force	normal reaction of a wire on a bead, tension in a string
Statistical physics	conjugate intensive quantity	\(\beta = 1/k_BT\) (energy constraint)
Chemistry	chemical potential \(\mu\)	particle-number constraint
Economics	shadow price	budget constraint = marginal utility

The mechanics row is the most intuitive: for a bead constrained to a slope, the \(\lambda\) solved from the Lagrange equations is precisely the normal reaction force the slope exerts on the bead—“the reaction force is perpendicular to the contact surface” and “\(\nabla g\) is perpendicular to the constraint set” are the same statement. The thermodynamic table of intensive–extensive pairs (\(T\)–\(S\), \(p\)–\(V\), \(\mu\)–\(N\)) is essentially a registry of Lagrange multipliers: each intensive quantity is the multiplier of some conserved-quantity constraint, which is also why they are inherently size-independent—a multiplier is a ratio \(\partial f^*/\partial c\). The “Boltzmann distribution” chapter will show the first row: temperature is the shadow price of the energy constraint, \(\partial S/\partial E = 1/T\).

Generalizing to distributions: in the maximum-entropy problem the “variables” are not two numbers but the density value \(p(x)\) at every point—infinitely many variables, one per \(x\). The recipe is unchanged (this step is called the calculus of variations), and the next subsection does exactly this. This method is everywhere in physics, because physical problems are almost always “some quantity is extremized + conservation laws as constraints,” and the multiplier is often a named physical quantity: temperature (energy constraint), chemical potential (particle-number constraint), pressure (volume constraint).

2.6 Lagrange multiplier derivation: the general form of the distribution

Maximize \(H = -\int p \ln p\, dx\) subject to the constraints:

1. Normalization: \(\int p(x) dx = 1\)
2. The expectation of some function is fixed: \(\int p(x) f(x) dx = \langle f \rangle\)

(Constraint 2 is not derived—it is the input to the problem: you measured, or a conservation law fixed, some average quantity, and “knowing an average” written mathematically is exactly an expectation—\(f(x)\) specifies which quantity is measured, and \(\langle f\rangle\) is the number measured. Compare:

• measured the sample mean \(\bar x\): take \(f(x) = x\), constraint \(\int p(x)\,x\,dx = \bar x\), i.e. \(\langle f\rangle = \bar x\);
• measured the variance \(s^2\) (together with the mean \(\bar x\)): add \(f(x) = x^2\), constraint \(\int p(x)\,x^2\,dx = s^2 + \bar x^2\) (second moment), i.e. \(\langle f\rangle = s^2 + \bar x^2\);
• energy conserved, each particle getting \(\bar E\) on average: take \(f(x) = E(x)\), constraint \(\int p(x)\,E(x)\,dx = \bar E\), i.e. \(\langle f\rangle = \bar E\).

Both \(f\) and \(\langle f\rangle\) are determined by known information; the maximum-entropy principle only ensures no extra assumptions beyond the known—so a different \(f\) means a different distribution.)

Construct the Lagrangian:

\[L = -\int p \ln p\, dx - \lambda_0\!\left(\int p\, dx - 1\right) - \lambda_1\!\left(\int p\, f(x)\, dx - \langle f\rangle\right)\]

Take the variation with respect to \(p(x)\) and set it to zero:

What does “take the variation” mean? Ordinary calculus optimizes a number; here we optimize an entire function (finding which \(p\) maximizes \(H[p]\))—this is the calculus of variations. In practice one rule suffices: differentiate the integrand \(F(p)\) with respect to \(p\) as if it were an ordinary variable,

\[\frac{\delta}{\delta p(x)} \int F(p)\, dx = \frac{\partial F}{\partial p}\]

Applied to the entropy term: \(F = p\ln p\), \(\frac{\partial}{\partial p}(p\ln p) = \ln p + 1\) (product rule, the \(+1\) coming from \(p\cdot\frac{1}{p}\)). So \(-\int p\ln p\) contributes \(-\ln p - 1\), the normalization term contributes \(-\lambda_0\), and the expectation term contributes \(-\lambda_1 f(x)\); adding them gives the equation below. (The \(+1\) is finally absorbed into \(Z = e^{1+\lambda_0}\) and does not affect the solution.)

\[\frac{\delta L}{\delta p} = -\ln p(x) - 1 - \lambda_0 - \lambda_1 f(x) = 0\]

Solving:

\[\boxed{p(x) = \frac{1}{Z}\, e^{-\lambda_1 f(x)}}\]

where \(Z = e^{1+\lambda_0}\) is the normalization constant (partition function), and \(\lambda_1\) is determined by the constraint \(\langle f\rangle\).

That is the entire derivation. Every distribution comes out of this one formula; the only difference is what \(f(x)\) is.

2.7 Constraints determine the distribution

par(mfrow = c(2, 2), mar = c(4, 4, 3.5, 1))

# ── 1. No constraint (bounded interval) → Uniform ──
x1 <- seq(-0.5, 1.5, length.out = 300)
plot(x1, dunif(x1, 0, 1), type = "l", col = "steelblue", lwd = 2.5,
     main = "No constraint (bounded)\n=> Uniform",
     xlab = "x", ylab = "p(x)", ylim = c(0, 1.5))
text(0.5, 1.3, "f(x) = none", font = 3, col = "gray40")
text(0.5, 0.4, "Maximum disorder\non [0, 1]", col = "gray40", cex = 0.95)

# ── 2. Fixed mean (nonnegative) → Exponential ──
x2 <- seq(0, 8, length.out = 300)
plot(x2, dexp(x2, rate = 1), type = "l", col = "steelblue", lwd = 2.5,
     main = "Constraint: E[X] = mu, X >= 0\n=> Exponential(1/mu)",
     xlab = "x", ylab = "p(x)", ylim = c(0, 1.1))
text(4, 0.8, expression(f(x) == x), font = 3, col = "gray40", cex = 1.1)
text(4, 0.6, expression(lambda[1] == 1/mu), col = "tomato", cex = 0.95)

# ── 3. Fixed mean + variance → Normal ──
x3 <- seq(-4, 4, length.out = 300)
plot(x3, dnorm(x3), type = "l", col = "steelblue", lwd = 2.5,
     main = "Constraint: E[X]=mu, Var(X)=sigma^2\n=> Normal(mu, sigma^2)",
     xlab = "x", ylab = "p(x)", ylim = c(0, 0.45))
text(2, 0.38, expression(f(x) == x^2), font = 3, col = "gray40", cex = 1.1)
text(2, 0.30, expression(lambda[2] == 1/(2*sigma^2)), col = "tomato", cex = 0.95)

# ── 4. Fixed log-mean → Power law / Pareto ──
x4 <- seq(1, 8, length.out = 300)
lambda4 <- 2.5
xmin4 <- 1
p4 <- (lambda4 - 1) * xmin4^(lambda4 - 1) * x4^(-lambda4)
plot(x4, p4, type = "l", col = "steelblue", lwd = 2.5,
     main = "Constraint: E[log X] fixed, X >= xmin\n=> Power law / Pareto",
     xlab = "x", ylab = "p(x)", ylim = c(0, 1.6))
text(4.8, 1.20, expression(f(x) == log(x)), font = 3, col = "gray40", cex = 1.1)
text(4.8, 0.95, expression(lambda > 1), col = "tomato", cex = 0.95)
text(4.8, 0.70, expression(p(x) %prop% x^-lambda), col = "gray40", cex = 0.95)

Constraint	\(f(x)\)	\(\lambda\)	Maximum-entropy distribution \(p(x)\)
No constraint (bounded)	—	—	Uniform: \(\dfrac{1}{b-a}\)
\(E[X] = \mu\), \(X \geq 0\)	\(x\)	\(1/\mu\)	Exponential: \(\lambda e^{-\lambda x}\)
\(E[\ln X]\) fixed, \(X \geq x_{\min}\)	\(\ln x\)	\(\lambda > 1\)	Power law: \((\lambda-1)x_{\min}^{\lambda-1}x^{-\lambda}\)
\(E[X] = \mu\), \(\text{Var}(X) = \sigma^2\)	\(x^2\)	\(1/(2\sigma^2)\)	Normal: \(\dfrac{1}{\sqrt{2\pi}\,\sigma}\,e^{-(x-\mu)^2/(2\sigma^2)}\)
\(E[\text{Energy}] = \langle E \rangle\)	\(E\)	\(\beta = 1/k_BT\)	Boltzmann: \(\beta\, e^{-\beta E}\)

All these distributions are different solutions of the same equation. The only difference is what you “know” (the constraints). Here “constraint” includes not only expectation conditions like mean, variance, and energy, but also the support itself: the uniform distribution first requires knowing the variable is confined to \([a,b]\), the exponential first requires \(X\ge0\), and the power law first requires \(X\ge x_{\min}>0\). The maximum-entropy problem is never “find a distribution out of thin air,” but rather:

\[\text{given support + given some statistics} \quad\Longrightarrow\quad \text{maximum entropy within that information}\]

---

3 Exponential distribution = maximum entropy + nonnegativity + mean constraint

3.1 Derivation

Constraints: \(X \geq 0\), \(E[X] = \mu\). Here \(X\ge0\) is itself a support constraint: waiting times, lifetimes, heights, energies cannot be negative; fixing the mean on this half-line then gives the exponential distribution. Plugging into the general solution:

\[p(x) = \frac{1}{Z} e^{-\lambda x} = \lambda e^{-\lambda x}, \quad \lambda = 1/\mu\]

Normalization constant \(Z = 1/\lambda\). Differential entropy \(H = -\int_0^\infty p\ln p\,dx\):

\[ \begin{aligned} H &= -\int_0^\infty \lambda e^{-\lambda x}\,\ln\!\bigl(\lambda e^{-\lambda x}\bigr)\,dx \\[2pt] &= -\int_0^\infty \lambda e^{-\lambda x}\,(\ln\lambda - \lambda x)\,dx && \ln\bigl(\lambda e^{-\lambda x}\bigr) = \ln\lambda - \lambda x \\[2pt] &= -\ln\lambda \int_0^\infty p\,dx \;+\; \lambda \int_0^\infty x\,p\,dx && \text{split up, }p = \lambda e^{-\lambda x} \\[2pt] &= -\ln\lambda \cdot 1 \;+\; \lambda \cdot \tfrac{1}{\lambda} && \textstyle\int p\,dx = 1,\ \int x\,p\,dx = E[X] = 1/\lambda \\[2pt] &= 1 - \ln\lambda = 1 + \ln\mu && \mu = 1/\lambda \end{aligned} \]

The unit is nats (using \(\ln\)). The larger the mean, the more spread out the distribution, and \(H\) grows accordingly; when \(\mu < 1/e\), \(H < 0\)—differential entropy can be negative, unlike discrete entropy.

3.2 Numerical check: among the candidates, the exponential has the highest entropy

Compare entropies among all nonnegative distributions with mean \(= 2\):

# Numerically compute differential entropy on a [0, max] grid
x_pos <- seq(0.001, 20, length.out = 8000)
dx_p  <- diff(x_pos)[1]

h_pos <- function(p_vals) {
  p_vals <- pmax(p_vals, 1e-300)
  p_norm <- p_vals / (sum(p_vals) * dx_p)
  -sum(p_norm * log(p_norm)) * dx_p
}

mu <- 2   # fixed mean

# Various distributions with mean = 2
dists <- list(
  "Exp(1/2)"         = dexp(x_pos, rate = 1/mu),
  "Gamma(2, 1)"      = dgamma(x_pos, shape = 2,  rate = 2/mu),
  "Gamma(5, 2.5)"    = dgamma(x_pos, shape = 5,  rate = 5/mu),
  "Gamma(20, 10)"    = dgamma(x_pos, shape = 20, rate = 20/mu),
  "Weibull(k=2)"     = dweibull(x_pos, shape = 2,
                                scale = mu / gamma(1 + 1/2)),
  "LogNormal"        = dlnorm(x_pos, meanlog = log(mu) - 0.5, sdlog = 1),
  "Uniform[0,4]"     = dunif(x_pos, 0, 2 * mu)   # uniform with mean 2: entropy is actually smaller
)

entropies <- sapply(dists, h_pos)
h_theory  <- 1 + log(mu)   # analytic entropy of Exp(1/mu)

par(mfrow = c(1, 2), mar = c(4, 7, 3.5, 1))

# Left: density comparison
plot(NA, xlim = c(0, 10), ylim = c(0, 0.55),
     main = "Distributions with same mean = 2",
     xlab = "x", ylab = "Density")
dist_cols <- c("tomato", "steelblue", "forestgreen", "orange", "purple", "brown", "gray40")
for (i in seq_along(dists)) {
  lines(x_pos, dists[[i]], col = dist_cols[i], lwd = 2,
        lty = ifelse(i == 1, 1, 2))
}
legend("topright", names(dists), col = dist_cols, lwd = 2,
       lty = c(1, rep(2, 6)), bty = "n", cex = 0.95)

# Right: entropy comparison (horizontal bar chart)
barplot(rev(entropies), horiz = TRUE, col = rev(dist_cols),
        border = "white", las = 1, xlim = c(0, max(entropies) * 1.28),
        main = "Entropy comparison (mean = 2)",
        xlab = "Differential entropy H")
abline(v = h_theory, col = "tomato", lwd = 2, lty = 2)
text(h_theory * 1.045, 4.3, paste("Exp theoretical =", round(h_theory, 3)),
     col = "tomato", adj = 0.5, srt = 90, cex = 0.95)

The red bar for Exp(1/2) is longest—among all nonnegative distributions with mean \(= 2\), the exponential really does have the maximum entropy.

As \(k\) grows (Gamma(k) becoming more normal in shape), its entropy falls, because the distribution becomes more “certain.”

Note the bottom bar, Uniform[0,4] (also mean 2): its entropy \(\ln 4 \approx 1.386\) is smaller than the exponential’s \(1.693\). This looks counterintuitive—isn’t “uniform = maximum entropy” the usual saying? That statement holds only when the support is fixed and bounded, with no other constraint (“\(X\) must lie in \([0,4]\), nothing else known” → uniform). Here the support is all of \([0,\infty)\) and the constraint is the mean: \(U[0,4]\) forces the probability to be hard zero for \(x>4\), which amounts to secretly asserting “\(X\) never exceeds 4”—extra information you do not actually have, which instead costs entropy. The exponential sets no upper bound, fills the whole half-line, and assumes nothing beyond the mean, so its entropy is higher.

3.3 Rigorous proof: at fixed mean, the exponential maximizes entropy (not just “tallest in the figure”)

The figure compared only 7 candidates. The proposition to prove is: among all nonnegative distributions with mean fixed at \(\mu\), the exponential has maximum entropy (“fixed mean” cannot be dropped—without a constraint the entropy has no upper bound). Use the Gibbs inequality: for any density \(p\) on \([0,\infty)\) with mean \(\mu\), take the exponential \(q(x)=\tfrac1\mu e^{-x/\mu}\) with the same mean; its KL divergence is nonnegative (\(\int p\ln\frac{p}{q}\ge 0\)), which expands to:

(This is “guess the answer, then verify”: \(q\) is not chosen out of thin air—it is exactly the candidate from the variational derivation in 2.6; from another angle, it is also the unique \(q\) that collapses the cross-entropy term below into a constant, see “Details.”)

\[ \begin{aligned} 0 \le \int p\ln\frac{p}{q}\,dx &= \int p\ln p\,dx - \int p\ln q\,dx && \ln\tfrac{p}{q} = \ln p - \ln q \\[2pt] &= -H(p) - \int p\ln q\,dx && H(p) = -\!\int p\ln p\,dx \\[2pt] &= -H(p) - \int p\Bigl(-\ln\mu - \tfrac{x}{\mu}\Bigr)dx && \ln q = -\ln\mu - x/\mu \\[2pt] &= -H(p) + \ln\mu + \tfrac{1}{\mu}\!\int x\,p\,dx && \textstyle\int x\,p\,dx = \mu \\[2pt] &= -H(p) + \ln\mu + 1 \end{aligned} \]

Rearranging gives \(H(p) \le 1 + \ln\mu = H(q)\), with equality only when \(p = q\). At fixed mean the exponential is the unique maximum-entropy solution—not merely the tallest by coincidence in the figure.

The division of labor between the two roads (a key point throughout this article): finding the maximum-entropy distribution takes two steps—the variation/Lagrange method of 2.6 “finds” the candidate \(p\propto e^{-\lambda f}\) (but variation only gives a stationary point, not automatically a guarantee of a global maximum); the Gibbs inequality here “proves” it is the unique global maximum—for all distributions with the same constraint \(H(p)\le H(q)\), with equality only when \(p=q\), so any other distribution is strictly smaller and the exponential is the unique maximizer (the price is needing the candidate as input). In a phrase: variation finds it, Gibbs certifies it.

The “Details” below explain the KL divergence used here and why the exponential is chosen as \(q\). If you only want the conclusion, skip straight to the next section, “Memorylessness.”

3.3.1 Details: KL divergence and two choices

What KL divergence is. The KL divergence between two densities \(p\), \(q\) is \(D(p\|q) = \int p\ln\frac{p}{q}\,dx\), read as “the average extra log-cost of describing the true distribution \(p\) using \(q\).” This “cost” has a concrete backdrop: Shannon coding assigns \(-\log_2 p(x)\) bits to a symbol of probability \(p(x)\), and coding by the true \(p\) gives the shortest average code length (\(=H(p)\)); if you mistakenly believe the distribution is \(q\) and build a code from \(q\), the extra bits spent per symbol on average is exactly \(D(p\|q)\)—using the wrong model is never cheaper, so it is \(\ge 0\). The cross-entropy loss in machine learning is precisely this quantity: training a classifier = making the model distribution \(q\) approach the true labels \(p\), i.e. minimizing \(D(p\|q)\). It is always nonnegative, with equality only when \(p=q\):

Watch next (optional). But what is cross-entropy? | Compression is Intelligence Part 2 (3Blue1Brown) connects cross-entropy, coding cost, and machine-learning loss; it is most useful after the coding interpretation above.

\[-D(p\|q) = \int p\ln\frac{q}{p}\,dx = E_p\Bigl[\ln\tfrac{q}{p}\Bigr] \;\overset{\text{Jensen}}{\le}\; \ln E_p\Bigl[\tfrac{q}{p}\Bigr] = \ln\!\int p\cdot\frac{q}{p}\,dx = \ln\!\int q\,dx = 0\]

That Jensen step says: for a concave function (∩-shaped, opening downward, like \(\ln\)), “average first, then apply the function” \(\ge\) “apply the function first, then average”, i.e. \(\varphi(E[X]) \ge E[\varphi(X)]\). Geometrically, any chord of a concave function lies below the curve, so the curve’s height at the “average point” exceeds the average of the endpoint values. (Note the terminology: a concave function is ∩-shaped, while a convex function like \(x^2\) is ∪-shaped. Watch the examples, not the words.) \(\ln\) is concave everywhere (\(\ln'' x = -1/x^2 < 0\)), so \(\ln E[Y] \ge E[\ln Y]\); here take \(Y = q/p\) with expectation under \(p\), and pulling \(\ln\) outside the expectation only increases the value, giving \(\le \ln E_p[q/p] = \ln\!\int q = \ln 1 = 0\). Hence \(-D \le 0\), i.e. \(D \ge 0\).

Why use KL. It is always \(\ge 0\), a ready-made tool for prying an “\(=\)” into a “\(\le\).” Splitting the first line above gives \(H(p) \le -\int p\ln q\), whose right side is the cross-entropy—entropy does not exceed the cross-entropy computed with any \(q\). The rest is just evaluating the cross-entropy.

Why choose the exponential. For a general \(q\) the cross-entropy is intractable. The exponential’s \(\ln q = -\ln\mu - x/\mu\) is linear in \(x\), so the cross-entropy involves only the mean \(\int x\,p = \mu\) (already fixed)—it no longer depends on the other details of \(p\) and collapses to the same constant \(\ln\mu+1 = H(q)\) for every \(p\) with mean \(\mu\). What does the work is exactly the linearity of \(f(x)=x\) (the \(f\) in the 2.6 general solution).

General recipe. Change the constraint, change the \(q\), and the same proof goes through: constraint \(\int f\,p = c\) → take \(q \propto e^{-\lambda f}\) (\(\ln q\) linear in \(f\)) → the cross-entropy sees only \(\int f\,p = c\) → collapses to a constant. \(f=x\to\) exponential, \(f=x^2\to\) normal, \(f=E\to\) Boltzmann. This shows the “variation finds it, Gibbs certifies it” division of labor holds for every distribution, not just the exponential.

3.4 Property: memorylessness

The exponential distribution has a further special property independent of maximum entropy: it is the unique continuous distribution satisfying memorylessness.

\[P(X > s + t \mid X > s) = P(X > t) \quad \forall\, s, t > 0\]

X_mem    <- rexp(N, rate = 1)
s        <- 1.5
residual <- X_mem[X_mem > s] - s

par(mfrow = c(1, 3), mar = c(4, 4, 3.5, 1))

hist(X_mem, breaks = 100, freq = FALSE, col = "steelblue", border = "white",
     main = "Original Exp(1)", xlab = "x", xlim = c(0, 8), ylim = c(0, 1.05))
curve(dexp(x, 1), add = TRUE, col = "tomato", lwd = 2)

hist(residual, breaks = 100, freq = FALSE, col = "steelblue", border = "white",
     main = paste0("Residual  (given X > ", s, ")"),
     xlab = "residual", xlim = c(0, 8), ylim = c(0, 1.05))
curve(dexp(x, 1), add = TRUE, col = "tomato", lwd = 2)
legend("topright", "Same Exp(1)!", col = "tomato", lwd = 2, bty = "n")

t_vals      <- seq(0, 5, by = 0.5)
p_original  <- pexp(t_vals, lower.tail = FALSE)
p_condition <- sapply(t_vals, function(t) mean(residual > t))
plot(t_vals, p_original, type = "l", col = "steelblue", lwd = 2,
     main = "Memoryless property verified",
     xlab = "t", ylab = "P(> t)")
points(t_vals, p_condition, col = "tomato", pch = 16, cex = 1.2)
legend("topright",
       c("P(X > t)  [unconditional]",
         paste0("P(X > t | X > ", s, ")")),
       col = c("steelblue", "tomato"), lwd = c(2, NA), pch = c(NA, 16), bty = "n")

Mathematical essence: \(P(X > t) = e^{-\lambda t}\)—the exponential function’s “for every small increase, the probability drops by a fixed proportion” directly implies memorylessness.

3.5 Physical examples: the isothermal atmosphere and radioactive decay

The barometric formula. An air molecule’s gravitational potential energy at height \(h\) is \(E = mgh\)—energy is linear in height. This is the physical embodiment of “linear constraint → exponential,” and the Boltzmann factor gives directly:

\[P(h) \propto e^{-mgh/k_BT}\]

Scale height \(H_{\text{scale}} = k_BT/mg \approx 8.5\) km (Earth’s atmosphere, \(T \approx 288\) K; atmospheric physics normally writes plain \(H\), subscripted here to avoid clashing with the \(H\) of entropy): every 8.5 km of altitude drops the pressure to \(1/e\). This also foreshadows the “Boltzmann distribution” chapter—the same formula, with energy levels there.

Radioactive decay. A nucleus does not age: a \(^{14}\)C nucleus that has existed for ten thousand years has exactly the same probability of decaying in the next second as a freshly created one—this is the memorylessness of the previous subsection, so lifetimes can only follow the exponential distribution. The half-life \(t_{1/2} = \ln 2/\lambda\) halves the surviving fraction at equal intervals: \(1/2, 1/4, 1/8, \dots\). First-order reaction kinetics in chemistry (\(-d[A]/dt = k[A]\)) and fluorescence-lifetime decay are the same mathematics.

## ── Physical examples: isothermal atmosphere + radioactive decay ──────────

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))

# Left: barometric formula P(h) = exp(-h/H), H = kT/mg ≈ 8.5 km
h_seq <- seq(0, 35, length.out = 400)
H_scale <- 8.5   # km, scale height of Earth's atmosphere
plot(h_seq, exp(-h_seq / H_scale), type = "l", lwd = 3, col = "steelblue",
     main = "Barometric formula:  P(h) = exp(-mgh / kT)",
     xlab = "Altitude h (km)", ylab = "Pressure relative to sea level",
     cex.lab = 1.35, cex.main = 1.1)
landmarks <- data.frame(h = c(8.85, 11, 30))
points(landmarks$h, exp(-landmarks$h / H_scale), pch = 16, col = "tomato", cex = 1.4)
text(9.3,  0.47, "Everest 8.8 km:  35%",      col = "tomato", adj = 0, cex = 1.05)
text(11.5, 0.20, "cruise 11 km:  27%",        col = "tomato", adj = 0, cex = 1.05)
text(29.5, 0.12, "stratosphere 30 km:  3%",   col = "tomato", adj = 1, cex = 1.05)
text(22, 0.75, "scale height\nH_scale = kT/mg = 8.5 km", col = "gray40", cex = 1.1)

# Right: radioactive decay, survival curve + equally spaced half-lives
lam    <- 1                       # decay constant
t_half <- log(2) / lam
lifetimes <- rexp(N, rate = lam)  # simulate lifetimes of N nuclei
t_seq  <- seq(0, 5, length.out = 300)
surv_sim <- sapply(t_seq, function(t) mean(lifetimes > t))
plot(t_seq, surv_sim, type = "l", lwd = 3, col = "steelblue",
     main = "Radioactive decay:  survival = exp(-t / tau)",
     xlab = "Time t", ylab = "Fraction of nuclei surviving",
     cex.lab = 1.35, cex.main = 1.1)
curve(exp(-lam * x), add = TRUE, col = "tomato", lwd = 2, lty = 2)
halves <- t_half * (1:4)
points(halves, 0.5^(1:4), pch = 16, col = "forestgreen", cex = 1.4)
segments(halves, 0, halves, 0.5^(1:4), col = "forestgreen", lty = 3)
text(halves, 0.5^(1:4) + 0.06,
     c("1/2", "1/4", "1/8", "1/16"), col = "forestgreen", cex = 1.15)
legend("topright",
       c("simulated 100k nuclei", "exp(-t/tau) theory",
         paste0("half-life = ln2/lambda = ", round(t_half, 2))),
       col = c("steelblue", "tomato", "forestgreen"),
       lwd = c(3, 2, NA), lty = c(1, 2, NA), pch = c(NA, NA, 16),
       bty = "n", cex = 1.0)

3.5.1 What decides the shape is the form of \(f(x)\)

Both of these examples yield the exponential distribution, and the root cause in both is the form of \(f(x)\) in the general solution \(p(x) \propto e^{-\lambda f(x)}\)—whichever degree \(\ln p\) is as a function of \(x\) decides which distribution you get:

\[f(x) = x\ (\text{linear}) \to \text{exponential}, \qquad f(x) = x^2 \to \text{half-normal}, \qquad \dots\]

“Constraining only the mean” written mathematically is \(\int x\,p\,dx = \mu\), i.e. \(f(x) = x\) linear—this step forces the exponential. Change the constrained quantity and the shape changes: on the same \([0,\infty)\), constraining the second moment (\(f(x) = x^2\)) yields the half-normal \(p \propto e^{-\lambda x^2}\). So what decides it is the functional form of the constrained quantity in \(x\), not the number of constraints or parameters—“having only one scale parameter” alone is not enough to pin down the exponential (Rayleigh and half-normal are also one-parameter distributions on \([0,\infty)\)).

The linearity in the two examples comes from different sources:

• Isothermal atmosphere: the potential energy \(E = mgh\) is linear in height → exponential. If the potential energy were \(\propto h^2\), the same temperature parameter would give a Gaussian.
• Radioactive decay: the corresponding linearity condition is a constant hazard rate (the per-unit-time decay probability does not change with age), i.e. memorylessness—\(\ln(\text{survival fraction}) = -\lambda t\) is linear in \(t\).

In a phrase: exponential \(\Leftrightarrow \ln p\) linear in \(x \Leftrightarrow\) constant hazard rate \(\Leftrightarrow\) memoryless. The difference among exponential, normal, and Boltzmann is just whether \(\ln p\) is a linear, quadratic, or energy function of \(x\)—back to the general solution of Section 2.6.

---

4 Power-law distribution = maximum entropy + log-mean constraint

The power-law distribution \(p(x) \propto x^{-\lambda}\) is everywhere in nature—earthquake energy, city populations, word frequencies, wealth, website link counts, species abundance… yet the classroom rarely explains why it appears, and even more rarely its intrinsic relation to the exponential and the normal. The answer is surprisingly simple: the same maximum-entropy framework spits out the power law once you change the constraint function. The exponential constrained the arithmetic mean \(\langle x\rangle\), the normal constrains the variance, and here the constraint is the log-mean \(\langle\ln x\rangle\), equivalent to fixing the geometric mean \(e^{E[\ln X]}\).

4.1 Derivation

Constraints: \(x \ge x_{\min} > 0\), and the mean of the logarithm is fixed:

\[E[\ln X] = m \qquad (m > \ln x_{\min})\]

Take the constraint function \(f(x) = \ln x\) and plug into the general solution:

\[p(x) = \frac{1}{Z}\, e^{-\lambda \ln x} = \frac{1}{Z}\, x^{-\lambda}\]

This is the power law. \(\ln p = -\lambda\ln x - \ln Z\) is linear in \(\ln x\)—so the power law is a straight line in log-log coordinates, with slope \(-\lambda\). This log-log straight line is the literal picture of “\(\ln p\) linear in the constrained quantity \(\ln x\),” exactly parallel to the exponential being a straight line on lin-log axes.

Normalization requires a lower bound \(x_{\min}\). Like the exponential, the power law is not a free distribution on the whole real line: the exponential’s natural support is \([0,\infty)\), and the power law’s natural support is \([x_{\min},\infty)\). The difference is that the exponential has no convergence problem near 0, whereas \(x^{-\lambda}\) diverges if extended all the way to 0, so the power law must have a strictly positive lower bound \(x_{\min}>0\). Then

\[Z = \int_{x_{\min}}^{\infty} x^{-\lambda}\,dx = \frac{x_{\min}^{\,1-\lambda}}{\lambda - 1} \quad(\text{requires }\lambda > 1\text{ to converge})\]

so the normalized power law (i.e. the Pareto distribution) is

\[p(x) = (\lambda-1)\, x_{\min}^{\,\lambda-1}\, x^{-\lambda}, \qquad x \ge x_{\min}\]

Notation convention: throughout this article \(\lambda\) is always the density exponent (\(p\propto x^{-\lambda}\), i.e. the multiplier in the maximum-entropy general solution). What the literature usually calls the tail exponent is the \(\lambda-1\) of the survival function \(P(X>x)\propto x^{-(\lambda-1)}\)—the two differ by 1. The “survival exponent \(\alpha\)” in the wealth example later is exactly \(\lambda-1\); don’t confuse them.

\(\lambda > 1\) is the threshold for the power law to “survive” (otherwise it cannot integrate to a finite probability), of the same “support/convergence” kind as the exponential’s \(x\ge0\) or the normal’s needing no truncation. The constraint \(m=E[\ln X]\) then determines \(\lambda\):

\[ E[\ln X] = \ln x_{\min} + \frac{1}{\lambda - 1} \quad\Longrightarrow\quad \lambda = 1 + \frac{1}{m - \ln x_{\min}} \]

This step is quite illuminating: if \(m\) is only slightly larger than \(\ln x_{\min}\), the scale hugs the lower bound, \(\lambda\) is large, and the tail is thin; if \(m\) is very large, the geometric mean is far from the lower bound, \(\lambda\) approaches 1, and the tail is extremely heavy.

The differential entropy can likewise be computed directly:

\[ \begin{aligned} H &= -\int_{x_{\min}}^\infty p(x)\ln p(x)\,dx \\[2pt] &= -\ln(\lambda-1) -(\lambda-1)\ln x_{\min} + \lambda E[\ln X] \\[2pt] &= \ln x_{\min} - \ln(\lambda-1) + 1 + \frac{1}{\lambda-1} \end{aligned} \]

Note that \(\lambda\) here is the power exponent of the density: \(p(x)\propto x^{-\lambda}\). Therefore a finite \(E[X]\) requires \(\lambda>2\), and finite variance requires \(\lambda>3\). \(\lambda>1\) only guarantees “the total probability can be normalized,” not that the mean and variance exist.

4.2 The three thresholds of the power exponent: 1, 2, 3

The easiest thing to misunderstand about power laws is conflating “can be normalized,” “finite mean,” and “finite variance” into one thing. They are actually three different thresholds:

\(\lambda\) range	Total probability	Mean \(E[X]\)	Variance \(\mathrm{Var}(X)\)	Intuition
\(\lambda \le 1\)	not finite	—	—	not a valid probability distribution
\(1 < \lambda \le 2\)	finite	infinite	infinite	extreme values so strong that even the theoretical mean fails to exist
\(2 < \lambda \le 3\)	finite	finite	infinite	a typical average level exists, but fluctuations are dominated by extreme values
\(\lambda > 3\)	finite	finite	finite	still heavy-tailed, but mean and fluctuation are describable by ordinary moments

Why these thresholds? Because the \(r\)-th moment requires

\[E[X^r] = \int_{x_{\min}}^\infty x^r p(x)\,dx \propto \int_{x_{\min}}^\infty x^{r-\lambda}\,dx\]

and \(\int^\infty x^a dx\) converges only when \(a<-1\), so

\[E[X^r] < \infty \quad\Longleftrightarrow\quad r-\lambda < -1 \quad\Longleftrightarrow\quad \lambda > r+1\]

Setting \(r=0\) gives \(\lambda>1\) (total probability), \(r=1\) gives \(\lambda>2\) (mean), and \(r=2\) gives \(\lambda>3\) (second moment/variance).

## ── Power-exponent thresholds: normalization, mean, variance ───────────────

lam_grid <- seq(1.05, 5, length.out = 600)
xmin_thr <- 1
mean_thr <- ifelse(lam_grid > 2, (lam_grid - 1) / (lam_grid - 2) * xmin_thr, NA)
var_thr <- ifelse(lam_grid > 3,
                  (lam_grid - 1) * xmin_thr^2 /
                    ((lam_grid - 3) * (lam_grid - 2)^2),
                  NA)

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))

plot(NA, xlim = c(1, 5), ylim = c(0, 10),
     main = "Power-law moments switch on at thresholds",
     xlab = expression(lambda), ylab = "Moment value (xmin = 1)",
     cex.lab = 1.35, cex.main = 1.05)
rect(1, 0, 2, 10, col = adjustcolor("gray70", 0.25), border = NA)
rect(2, 0, 3, 10, col = adjustcolor("tomato", 0.12), border = NA)
lines(lam_grid, pmin(mean_thr, 10), col = "steelblue", lwd = 3)
lines(lam_grid, pmin(var_thr, 10), col = "purple", lwd = 3)
abline(v = c(1, 2, 3), col = "gray40", lty = 2)
text(1.5, 8.6, "not enough\nfor mean", col = "gray35", cex = 0.95)
text(2.5, 8.6, "mean finite,\nvariance\ninfinite", col = "tomato", cex = 0.95)
text(3.8, 5.6, "mean +\nvariance\nfinite", col = "gray35", cex = 0.95)
legend("topright", c("E[X]", "Var(X)"),
       col = c("steelblue", "purple"), lwd = 3, bty = "n", cex = 0.95)

x_tail <- seq(1, 1000, length.out = 2000)
plot(NA, xlim = c(1, 1000), ylim = c(1e-9, 2),
     log = "xy",
     main = "Same lower bound, different tail thickness",
     xlab = "x (log scale)", ylab = "p(x) (log scale)",
     cex.lab = 1.35, cex.main = 1.05)
lam_show <- c(1.5, 2.5, 3.5)
cols_show <- c("gray40", "tomato", "steelblue")
for (i in seq_along(lam_show)) {
  lam <- lam_show[i]
  lines(x_tail, (lam - 1) * x_tail^(-lam), col = cols_show[i], lwd = 3)
}
legend("topright",
       c("lambda = 1.5: mean infinite",
         "lambda = 2.5: mean finite, var infinite",
         "lambda = 3.5: mean + var finite"),
       col = cols_show, lwd = 3, bty = "n", cex = 0.95)

The range \(2<\lambda\le3\) is often discussed separately, because it has a very clear statistical meaning: a typical scale exists, but the fluctuations have no finite theoretical scale. Here one can speak of an average city size, average wealth, average connectivity; but the variance does not exist, and sample fluctuations are long dominated by a few extreme values. As the data grow, a new extremely large observation can still significantly rewrite the variance estimate.

The other two ranges have their own meanings. \(1<\lambda\le2\) is heavier-tailed: the total probability is already finite, but the theoretical mean does not exist and the sample average is very unstable. \(\lambda>3\) is much gentler: both mean and variance are finite, though the tail is still thicker than the exponential’s. Real data also often have finite system size, measurement lower bounds, upper-bound truncation, and mixed mechanisms, so an empirically estimated \(\lambda\) can land in different ranges; the point is not to memorize a fixed range but to read off which moments exist for the range \(\lambda\) falls in.

4.3 Numerical check: at fixed \(E[\ln X]\), the power law has the highest entropy

Take \(x_{\min}=1\) and fix \(E[\ln X]=0.5\). The maximum-entropy solution should be a Pareto-type power law with \(\lambda = 1 + 1/0.5 = 3\). Below we compare several candidate distributions, all with support \([1,\infty)\) and numerically tuned to the same \(E[\ln X]\).

## ── Entropy comparison of candidate distributions with fixed E[log X] ──────

xmin_pl <- 1
m_log   <- 0.5
lambda_pl <- 1 + 1 / (m_log - log(xmin_pl))

safe_entropy <- function(dfun, lower = xmin_pl, upper = Inf) {
  integrate(function(x) {
    p <- dfun(x)
    out <- numeric(length(p))
    keep <- p > 0
    out[keep] <- -p[keep] * log(p[keep])
    out
  }, lower = lower, upper = upper, rel.tol = 1e-9)$value
}

elog_shift_exp <- function(rate) {
  integrate(function(x) log(x) * rate * exp(-rate * (x - xmin_pl)),
            lower = xmin_pl, upper = Inf, rel.tol = 1e-9)$value
}
rate_shift_exp <- uniroot(function(r) elog_shift_exp(r) - m_log,
                          c(0.02, 20))$root

elog_shift_gamma <- function(rate, shape = 2) {
  integrate(function(x) log(x) * dgamma(x - xmin_pl, shape = shape, rate = rate),
            lower = xmin_pl, upper = Inf, rel.tol = 1e-9)$value
}
rate_shift_gamma <- uniroot(function(r) elog_shift_gamma(r) - m_log,
                            c(0.02, 20))$root

elog_shift_weibull <- function(scale, shape = 2) {
  integrate(function(x) log(x) * dweibull(x - xmin_pl, shape = shape, scale = scale),
            lower = xmin_pl, upper = Inf, rel.tol = 1e-9)$value
}
scale_shift_weib <- uniroot(function(s) elog_shift_weibull(s) - m_log,
                            c(0.05, 30))$root

elog_uniform <- function(b) (b * log(b) - b + xmin_pl) / (b - xmin_pl)
b_unif_pl <- uniroot(function(b) elog_uniform(b) - m_log, c(1.01, 20))$root

dists_pl <- list(
  "Power law" = function(x) (lambda_pl - 1) * xmin_pl^(lambda_pl - 1) * x^(-lambda_pl),
  "Shifted Exp" = function(x) rate_shift_exp * exp(-rate_shift_exp * (x - xmin_pl)),
  "Shifted Gamma(k=2)" = function(x) dgamma(x - xmin_pl, shape = 2, rate = rate_shift_gamma),
  "Shifted Weibull(k=2)" = function(x) dweibull(x - xmin_pl, shape = 2, scale = scale_shift_weib),
  "Uniform[1,b]" = function(x) dunif(x, xmin_pl, b_unif_pl)
)

ent_pl <- c(
  "Power law" = log(xmin_pl) - log(lambda_pl - 1) + 1 + 1 / (lambda_pl - 1),
  "Shifted Exp" = safe_entropy(dists_pl[["Shifted Exp"]]),
  "Shifted Gamma(k=2)" = safe_entropy(dists_pl[["Shifted Gamma(k=2)"]]),
  "Shifted Weibull(k=2)" = safe_entropy(dists_pl[["Shifted Weibull(k=2)"]]),
  "Uniform[1,b]" = log(b_unif_pl - xmin_pl)
)

elog_checks <- sapply(dists_pl, function(dfun) {
  integrate(function(x) log(x) * dfun(x),
            lower = xmin_pl, upper = Inf, rel.tol = 1e-8)$value
})

x_plot_pl <- seq(1, 12, length.out = 2000)
dist_cols_pl <- c("tomato", "steelblue", "forestgreen", "purple", "gray40")

par(mfrow = c(1, 2), mar = c(4, 7.2, 3.5, 1))

plot(NA, xlim = c(1, 10), ylim = c(0, 2.2),
     main = "Same constraint: E[log X] = 0.5",
     xlab = "x", ylab = "Density")
for (i in seq_along(dists_pl)) {
  lines(x_plot_pl, dists_pl[[i]](x_plot_pl), col = dist_cols_pl[i], lwd = 2,
        lty = ifelse(i == 1, 1, 2))
}
legend("topright", names(dists_pl), col = dist_cols_pl, lwd = 2,
       lty = c(1, rep(2, 4)), bty = "n", cex = 0.95)

barplot(rev(ent_pl), horiz = TRUE, col = rev(dist_cols_pl),
        border = "white", las = 1, xlim = c(0, max(ent_pl) * 1.18),
        main = "Entropy comparison",
        xlab = "Differential entropy H")
abline(v = ent_pl["Power law"], col = "tomato", lwd = 2, lty = 2)
text(ent_pl["Power law"] * 1.045, 3.0,
     paste("Power law theoretical =", round(ent_pl["Power law"], 3)),
     col = "tomato", adj = 0.5, srt = 90, cex = 0.95)

Every candidate distribution above satisfies the same constraint; numerical verification:

distribution	E_log_X_numeric
Power law	0.500000
Shifted Exp	0.500000
Shifted Gamma(k=2)	0.499999
Shifted Weibull(k=2)	0.500000
Uniform[1,b]	0.499998

The red bar for the power law is longest—with support \([x_{\min},\infty)\) and \(E[\ln X]\) fixed, the power-law distribution really does have the highest entropy.

Uniform looks “flat,” but it forces the probability to hard zero for \(x>b\), secretly adding the extra information “will never exceed \(b\)”; the shifted exponential/gamma/Weibull press the tail thinner than the power law, also saying more than the constraint. The power law fixes only the geometric scale, not the arithmetic mean, so it allows extreme values far more often than the exponential does.

4.4 Rigorous proof: at fixed log-mean, the power law maximizes entropy

This section is a skippable rigorous proof (isomorphic to the Gibbs argument in the exponential section). If you only want the conclusion—“at fixed \(E[\ln X]\) the power law maximizes entropy”—skip straight to the next section, “Scale-free.”

Use the Gibbs inequality, entirely isomorphic to the exponential case. For any density \(p\) on \([x_{\min},\infty)\) satisfying \(E_p[\ln X]=m\), take the power-law candidate with the same constraint

\[q(x) = (\lambda-1)x_{\min}^{\lambda-1}x^{-\lambda}, \qquad \lambda = 1 + \frac{1}{m-\ln x_{\min}}\]

The KL divergence is nonnegative:

\[ \begin{aligned} 0 \le \int p\ln\frac{p}{q}\,dx &= -H(p) - \int p\ln q\,dx \\[2pt] &= -H(p) - \int p\Bigl[\ln(\lambda-1)+(\lambda-1)\ln x_{\min}-\lambda\ln x\Bigr]dx \\[2pt] &= -H(p) -\ln(\lambda-1)-(\lambda-1)\ln x_{\min}+\lambda m \\[2pt] &= -H(p) + \ln x_{\min} - \ln(\lambda-1) + 1 + \frac{1}{\lambda-1} \end{aligned} \]

The last line is exactly \(H(q)\), so \(H(p)\le H(q)\), with equality only when \(p=q\). This is the same “variation finds it, Gibbs certifies it” routine: because \(\ln q\) is linear in \(\ln x\), the cross-entropy depends only on the fixed \(E[\ln X]\) and no longer on the other details of \(p\).

4.5 Power law = exponential on a log scale

This section reduces “power law” to a single sentence: it is just the exponential distribution, only living on a log scale. First see it geometrically, then nail it with one change of variable.

Geometric entry point: straighten its coordinate system and you read off its constraint. Take the log of the general solution—\(\ln\) both sides of \(p\propto e^{-\lambda f(x)}\):

\[\ln p(x) = -\lambda\, f(x) + \text{const}\]

\(\ln p\) is a linear function of \(f(x)\). Plotting \(\ln p\) on the vertical axis against \(f(x)\) on the horizontal, any maximum-entropy distribution is a straight line with slope \(-\lambda\). Conversely, “which coordinate system straightens the distribution” = “which \(f\) it constrains” = “which maximum-entropy distribution it is”—log-log and lin-log are not casually chosen plotting styles but “plotting \(\ln p\) against the constrained quantity.”

## Each maximum-entropy distribution is straight only when "horizontal axis = its f(x)"
x <- seq(1, 30, length.out = 1500)
p_exp   <- 0.15 * exp(-0.15 * x)                 # f = x
p_power <- (2.5 - 1) * x^(-2.5)                   # f = ln x, xmin = 1
p_gauss <- dnorm(x, 0, 8)                         # f = x^2 (mu = 0)
col_e <- "steelblue"; col_p <- "tomato"; col_g <- "forestgreen"
yl <- c(1e-4, 3)
par(mfrow = c(1, 3), mar = c(4.4, 4.6, 3.6, 1))

# A: axis x -> exponential is straight
plot(x, p_exp, type = "l", lwd = 4, col = col_e, log = "y", ylim = yl,
     main = "axis: x   =>  Exponential straight",
     xlab = "x", ylab = "density (log)", cex.lab = 1.25, cex.main = 1.05)
lines(x, p_power, col = col_p, lwd = 1.6, lty = 2)
lines(x, p_gauss, col = col_g, lwd = 1.6, lty = 2)
legend("topright", c("Exp  (f=x)", "Power", "Gauss"),
       col = c(col_e, col_p, col_g), lwd = c(4, 1.6, 1.6), lty = c(1, 2, 2),
       bty = "n", cex = 0.95)

# B: axis ln x -> power law is straight
plot(x, p_power, type = "l", lwd = 4, col = col_p, log = "xy", ylim = yl,
     main = "axis: ln x   =>  Power law straight",
     xlab = "x (log)", ylab = "density (log)", cex.lab = 1.25, cex.main = 1.05)
lines(x, p_exp, col = col_e, lwd = 1.6, lty = 2)
lines(x, p_gauss, col = col_g, lwd = 1.6, lty = 2)
legend("bottomleft", c("Power  (f=ln x)", "Exp", "Gauss"),
       col = c(col_p, col_e, col_g), lwd = c(4, 1.6, 1.6), lty = c(1, 2, 2),
       bty = "n", cex = 0.95)

# C: axis x^2 -> Gaussian is straight
xx <- x^2
plot(xx, p_gauss, type = "l", lwd = 4, col = col_g, log = "y", ylim = yl,
     main = "axis: x^2   =>  Gaussian straight",
     xlab = expression(x^2), ylab = "density (log)", cex.lab = 1.25, cex.main = 1.05)
lines(xx, p_exp, col = col_e, lwd = 1.6, lty = 2)
lines(xx, p_power, col = col_p, lwd = 1.6, lty = 2)
legend("topright", c("Gauss  (f=x^2)", "Exp", "Power"),
       col = c(col_g, col_e, col_p), lwd = c(4, 1.6, 1.6), lty = c(1, 2, 2),
       bty = "n", cex = 0.95)

Three panels, the same set of curves, only the horizontal axis changed: each distribution is straight only in the panel where “horizontal axis = its own \(f\).” Read as a table:

Distribution	Coordinates that straighten it (vert. \(\ln p\), horiz. \(\cdots\))	Slope	Constraint \(f(x)\)
Exponential	\(x\) (lin-log)	\(-\lambda\)	\(x\)
Power law	\(\ln x\) (log-log)	\(-\lambda\)	\(\ln x\)
Normal	\(x^2\)	\(-\tfrac{1}{2\sigma^2}\)	\(x^2\)

The punchline: the power law’s log-log straight line is exactly its being an exponential on a log scale. The power law being straight on log-log means \(\ln p\) is linear in \(\ln x\)—and “\(\ln p\) linear in some variable” is precisely the signature of an exponential distribution in that variable (the exponential is straight on lin-log). So the power law’s log-log straight line is equivalent to “treating \(\ln x\) as the new variable, it is a lin-log straight line = an exponential distribution.” A change of variable nails it: let \(Y=\ln(X/x_{\min})\ge 0\) (\(x=x_{\min}e^{y}\), \(dx/dy=x\)),

\[p_Y(y)\;\propto\;\bigl(x_{\min}e^{y}\bigr)^{-\lambda}\cdot x_{\min}e^{y} = e^{-(\lambda-1)y}\]

\(Y=\ln(X/x_{\min})\) follows exactly the exponential distribution with rate \(\lambda-1\):

\[\text{power law}(x)\;\Longleftrightarrow\;\text{exponential}(\ln x)\]

“\(\ln x\) is to the power law as \(x\) is to the exponential.” This also incidentally explains the constraint formula from the derivation, \(E[\ln X]=\ln x_{\min}+\tfrac{1}{\lambda-1}\)—the \(\tfrac{1}{\lambda-1}\) on the right is exactly the mean of an exponential distribution with rate \(\lambda-1\). Constraining the geometric mean \(\langle\ln x\rangle\) is exactly constraining the mean of “that exponential in log-space,” identical to the exponential’s constraint \(\langle x\rangle\). So the power law is not “another kind of” exponential—it is the exponential, just at a different scale.

From this, the parallel of memorylessness also emerges:

	Exponential	Power law
Memoryless form	\(P(X>s{+}t\mid X>s)=P(X>t)\)	\(P(X>ct\mid X>t)=P\!\left(\tfrac{X}{x_{\min}}>c\right)\)
Invariant operation	adding an amount \(t\) (translation)	multiplying by a factor \(c\) (scaling)
Characteristic scale	yes, \(1/\lambda\)	none (scale-free)

The exponential is additively memoryless (translation-invariant on a linear scale); the power law is multiplicatively memoryless (translation-invariant on a log scale = scale-invariant). Take one logarithm and multiplicative becomes additive, power law becomes exponential—the same “memorylessness,” just living on different scales. The next section draws this “multiplicative memorylessness” and proves the power law is the unique scale-free distribution.

4.6 Scale-free: draw “multiplicative memorylessness” and prove it is unique

The previous section already showed the power law is multiplicatively memoryless—\(P(X>ct\mid X>t)=c^{-(\lambda-1)}\), depending only on the multiplier \(c\), not the starting point \(t\). Here we do two things the previous section did not: draw it, and prove the power law is the unique such distribution.

## ── Scale-free: tail ratios depend on multiplier, not starting scale ─────

t_grid <- seq(1, 100, length.out = 600)
lambda_sf <- 2.5
mu_sf <- 10
multipliers <- c(2, 10)
cols_sf <- c("tomato", "steelblue")

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))

plot(NA, xlim = range(t_grid), ylim = c(1e-4, 1),
     log = "y",
     main = "Pareto: P(X > c t | X > t)",
     xlab = "Starting threshold t", ylab = "Conditional tail ratio",
     cex.lab = 1.35, cex.main = 1.05)
for (i in seq_along(multipliers)) {
  c_i <- multipliers[i]
  lines(t_grid, rep(c_i^(-(lambda_sf - 1)), length(t_grid)),
        col = cols_sf[i], lwd = 3)
}
legend("topright",
       paste0("c = ", multipliers, ": ratio = ",
              round(multipliers^(-(lambda_sf - 1)), 3)),
       col = cols_sf, lwd = 3, bty = "n", cex = 0.95)
text(48, 0.035, "flat lines:\nonly multiplier matters",
     col = "gray35", cex = 1.0)

plot(NA, xlim = range(t_grid), ylim = c(1e-40, 1),
     log = "y",
     main = "Exponential: P(X > c t | X > t)",
     xlab = "Starting threshold t", ylab = "Conditional tail ratio",
     cex.lab = 1.35, cex.main = 1.05)
for (i in seq_along(multipliers)) {
  c_i <- multipliers[i]
  lines(t_grid, exp(-(c_i - 1) * t_grid / mu_sf),
        col = cols_sf[i], lwd = 3)
}
legend("topright", paste0("c = ", multipliers),
       col = cols_sf, lwd = 3, bty = "n", cex = 0.95)
text(54, 1e-10, "curves fall with t:\nthere is a scale mu",
     col = "gray35", cex = 1.0)

The horizontal lines in the left panel are scale-freeness itself: for the power law, the probability of “doubling again” does not change with the starting point \(t\); for the exponential in the right panel, doubling again gets harder the further out you go (\(P(X>ct\mid X>t)=e^{-(c-1)t/\mu}\) depends on \(t\))—this is exactly the root of the exponential having a characteristic scale \(\mu\) and the power law having none.

Uniqueness: if a sufficiently regular positive function satisfies \(f(cx)=a(c)f(x)\) for any scaling, it can only be \(Kx^{-\lambda}\). So “exact scale-freeness” alone pins the distribution down to a power law—scale-free is not a property of the power law but another definition of it. (Empirical data satisfy it only approximately in the tail or over a finite range; truncated power laws and log-normals also look nearly straight on a finite log-log plot, so a straight line is a clue, not a proof—for rigorous determination see “What counts as a ‘real’ power law.”)

4.7 Scale-free = fractal = critical: three faces of the same coin

The just-proved “exact scale-freeness can only be a power law” steps straight into fractal geometry—because “scale-free” is exactly the definition of a fractal. This section joins the power law, fractals, and critical dynamics into one thing, and is also the chapter’s exit toward complex systems.

Self-similar ⟺ scale-invariant ⟺ power law. The defining property of a fractal is self-similarity: it looks the same magnified by any factor. “Shape unchanged after scaling” is \(f(cx)=a(c)f(x)\)—exactly the premise of the uniqueness theorem above, whose conclusion is that it can only be \(x^{-\lambda}\). So these three statements are three phrasings of the same proposition: fractal self-similarity, scale invariance, power law.

The fractal dimension is just a power exponent. How do you measure a fractal? Cover it with boxes of side \(\varepsilon\) and count the number of boxes \(N(\varepsilon)\); a fractal satisfies \(N(\varepsilon)\propto\varepsilon^{-D}\)—a power law, whose exponent \(D\) is the fractal dimension. Measuring a fractal = measuring a slope on log-log (the same routine as the previous section):

## Koch curve + box counting: the fractal dimension D is the slope of a log-log power law
koch_points <- function(level) {
  pts <- rbind(c(0,0), c(1,0))
  ang <- pi/3
  R <- matrix(c(cos(ang), sin(ang), -sin(ang), cos(ang)), 2, 2)
  for (it in seq_len(level)) {
    new <- pts[1, , drop=FALSE]
    for (i in 1:(nrow(pts)-1)) {
      A <- pts[i,]; B <- pts[i+1,]
      d <- (B - A)/3
      P1 <- A + d; P3 <- A + 2*d
      peak <- P1 + as.numeric(R %*% d)
      new <- rbind(new, P1, peak, P3, B)
    }
    pts <- new
  }
  pts
}
P <- koch_points(6)
dens <- do.call(rbind, lapply(1:(nrow(P)-1), function(i){
  A<-P[i,]; B<-P[i+1,]; t<-seq(0,1,length.out=6)
  cbind(A[1]+(B[1]-A[1])*t, A[2]+(B[2]-A[2])*t)
}))
boxcount <- function(pts, eps) nrow(unique(cbind(floor(pts[,1]/eps), floor(pts[,2]/eps))))
eps <- 1/2^(2:8)
Nbox <- sapply(eps, function(e) boxcount(dens, e))
D_box <- coef(lm(log(Nbox) ~ log(1/eps)))[2]

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.6, 1))
plot(P, type = "l", asp = 1, col = "steelblue", lwd = 1.5, axes = FALSE,
     xlab = "", ylab = "", main = "Koch curve: self-similar at every zoom")
plot(log10(1/eps), log10(Nbox), pch = 16, col = "tomato", cex = 1.5,
     main = "Box count  N(eps) ~ eps^(-D)",
     xlab = expression(log[10](1/epsilon)), ylab = expression(log[10]~N),
     cex.lab = 1.3, cex.main = 1.1)
abline(lm(log10(Nbox) ~ log10(1/eps)), col = "steelblue", lwd = 2.5)
legend("topleft",
       c(sprintf("slope = D = %.3f", D_box), sprintf("theory ln4/ln3 = %.3f", log(4)/log(3))),
       bty = "n", cex = 1.05, text.col = c("steelblue","gray30"))
text(mean(log10(1/eps)), min(log10(Nbox))+0.4,
     "fractal dimension\n= power-law exponent", col = "gray30", cex = 1.05)

The slope of the Koch-curve box count gives \(D=\ln4/\ln3\approx1.26\)—a non-integer dimension, precisely the incarnation of a power exponent. The Cantor set \(\ln2/\ln3\approx0.63\), Sierpinski \(\ln3/\ln2\approx1.58\), and coastlines \(L(\varepsilon)\propto\varepsilon^{1-D}\) (the finer you measure, the longer they get) are the same. So the power law’s exponent \(\lambda\) and the fractal dimension \(D\) are the same kind of thing—the exponent of “how some quantity scales with scale.”

Dynamical origins: chaos, SOC, criticality. Power laws/fractals are everywhere in nature because several classes of dynamics naturally converge to scale invariance: deterministic chaos settles onto strange attractors of fractal dimension (Lorenz, Hénon); self-organized criticality (the Bak–Tang–Wiesenfeld sandpile) makes a dissipative system automatically halt at the critical point, with avalanche sizes following a power law, no tuning needed—earthquakes (Gutenberg–Richter, the very row in the earlier table), forest fires, neuronal avalanches all belong here; at a phase-transition critical point the correlation length diverges, there is no characteristic scale, and correlations decay as a power law, with the renormalization group (RG) viewing the system as a fixed point of the scaling transformation—the deepest theory of “scale invariance \(\Rightarrow\) power law.”

The weld (the most beautiful part). This chapter gave the power law a maximum-entropy origin: “knowing nothing” under the scale-invariant measure \(dx/x\) (see the extension chapter later). And \(dx/x\) is exactly the unique measure invariant under \(x\to cx\)—it is the “scale invariance” of the RG/fractal language. Hence

\[\underbrace{\text{MaxEnt: maximally ignorant on a log scale}}_{\text{inference language}}\;=\;\underbrace{\text{RG: fixed point of the scaling transformation}}_{\text{dynamics language}}\]

Two languages, the same “no characteristic scale.” The uniqueness theorem above is its static shadow; the RG is its dynamical version.

A clue for biology. Fractals/scale-freeness are structural in biology: blood vessels, bronchi, and neuronal trees are all fractal branching networks; Kleiber’s metabolic law \(\propto M^{3/4}\) was derived by West–Brown–Enquist from fractal transport networks; the power law of neuronal avalanches supports the “critical brain” hypothesis. These are useful mechanistic clues, not a theorem that every living system sits at a critical point. (West et al. 1997; Beggs and Plenz 2003)

4.8 Combining: constraining \(\langle x\rangle\) and \(\langle\ln x\rangle\) at once → Gamma

The exponential fixes only the arithmetic mean \(\langle x\rangle\), the power law only the geometric mean \(\langle\ln x\rangle\). What if you fix both at once? Two constraints get two multipliers \(\lambda_1\) (paired with \(x\)) and \(\lambda_2\) (paired with \(\ln x\)); plug into the general solution:

\[p(x)\propto e^{-\lambda_1 x - \lambda_2\ln x} = x^{-\lambda_2}\,e^{-\lambda_1 x}\]

Power-law core \(\times\) exponential tail. On \((0,\infty)\) this is exactly the Gamma distribution (shape \(k=1-\lambda_2\), rate \(\lambda_1\)); on \([x_{\min},\infty)\) it is the literature’s “power law with exponential cutoff.” Two knobs, two limits:

• \(\lambda_1\to 0\) (release the arithmetic mean) → the exponential tail vanishes → pure power law \(x^{-\lambda_2}\);
• \(\lambda_2\to 0\) (release the geometric mean) → the power-law core vanishes → pure exponential \(e^{-\lambda_1 x}\).

So the exponential and the power law are not two species but two ends of one family: \(\lambda_2\) governs the small-to-medium-scale power-law core (log-log slope \(-\lambda_2\)), and \(\lambda_1\) governs the large-scale exponential cutoff (bending the tail down around \(x\sim 1/\lambda_1\)). In the next two sections, the subcritical branching of cell lineages and the Kesten steady state of wealth both have tails of exactly this “core + cutoff” form—they do not coincidentally “resemble” a Gamma, they are the maximum-entropy solution under two constraints. This is also the first appearance of “two constraints”: the next chapter, the normal, swaps the constraints for \(x\) and \(x^2\), again two multipliers.

## Gamma bridge: p ∝ x^{-λ2} e^{-λ1 x}. Two knobs: λ1 = exponential-tail cutoff, λ2 = power-law-core slope
xg <- exp(seq(log(1), log(3e4), length.out = 1500))
dens <- function(x, l2, l1) x^(-l2) * exp(-l1 * x)
norm_plot <- function(x, l2, l1) { y <- dens(x, l2, l1); y / max(y) }

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.6, 1))

# Left: fix core λ2=2, vary cutoff λ1
l1s <- c(0, 0.003, 0.05); colsL <- c("tomato", "purple", "steelblue")
plot(NA, xlim = c(1, 3e4), ylim = c(1e-12, 1.5), log = "xy",
     main = expression("Fix core "*lambda[2]*"=2, vary cutoff "*lambda[1]),
     xlab = "x (log)", ylab = "p(x) / peak (log)", cex.lab = 1.25, cex.main = 1.1)
for (i in seq_along(l1s)) lines(xg, norm_plot(xg, 2, l1s[i]), col = colsL[i], lwd = 3)
legend("topright",
       c(expression(lambda[1]*"=0  (pure power law)"),
         expression(lambda[1]*"=0.003"),
         expression(lambda[1]*"=0.05  (early cutoff)")),
       col = colsL, lwd = 3, bty = "n", cex = 0.95)
text(6, 3e-6, "straight core:\nslope -lambda2", col = "gray30", cex = 0.95)
text(6000, 3e-6, "exp cutoff\nknee ~ 1/lambda1", col = "gray30", cex = 0.95)

# Right: fix cutoff λ1=0.01, vary core λ2; λ2=0 = pure exponential
l2s <- c(0, 1, 2.5); colsR <- c("forestgreen", "purple", "tomato")
plot(NA, xlim = c(1, 3e4), ylim = c(1e-12, 1.5), log = "xy",
     main = expression("Fix cutoff "*lambda[1]*"=0.01, vary core "*lambda[2]),
     xlab = "x (log)", ylab = "p(x) / peak (log)", cex.lab = 1.25, cex.main = 1.1)
for (i in seq_along(l2s)) lines(xg, norm_plot(xg, l2s[i], 0.01), col = colsR[i], lwd = 3)
legend("topright",
       c(expression(lambda[2]*"=0  (pure exponential)"),
         expression(lambda[2]*"=1"),
         expression(lambda[2]*"=2.5  (steep core)")),
       col = colsR, lwd = 3, bty = "n", cex = 0.95)

Left: fixing the core \(\lambda_2\) and increasing \(\lambda_1\), the “knee” of the exponential cutoff moves from right to left, and at \(\lambda_1=0\) it reverts to a straight pure power law. Right: fixing the cutoff \(\lambda_1\) and increasing the core \(\lambda_2\), the core on log-log gets steeper, and at \(\lambda_2=0\) it reverts to a pure exponential. In reality, “finite-size truncated power laws” almost all live inside this figure.

4.9 Where real-world power exponents land

The preceding was all framework; here it lands. Some repeatedly studied power laws are listed by tail exponent \(\alpha\) (survival function \(P(X>x)\propto x^{-\alpha}\), i.e. density exponent \(\lambda-1\))—which moment range it falls in (see “The three thresholds of the power exponent” earlier) directly decides “whether the mean and variance still mean anything”:

Phenomenon	Quantity	Tail exponent \(\alpha\ (\approx)\)	Density \(\lambda=\alpha{+}1\)	Moment status
Zipf word frequency / city population	count, population	\(\sim 1\)	\(\sim 2\)	mean marginal, variance infinite
Wealth / income upper tail (Pareto)	wealth	\(1\text{–}2\)	\(2\text{–}3\)	often “mean finite, variance infinite”
Network degree distribution (scale-free nets)	node degree	\(2\text{–}3\)	\(3\text{–}4\)	mean finite, variance often infinite
Paper citations / website inlinks	citations, inlinks	\(2\text{–}3\)	\(3\text{–}4\)	as above
Earthquakes (Gutenberg–Richter)	released energy	\(\sim 2/3\)	\(\sim 5/3\)	even the energy mean fails to converge, extreme events dominate
Gene family size (genome evolution)	homolog-family members	\(\sim 2\)	\(\sim 3\)	gene duplication + divergence (preferential growth)

A one-line reminder: these exponents are all empirical estimates, approximate over a range. Rigorous statistical tests show that much of what is called a “power law” is not a clean power law over its whole range but a truncated power law, log-normal, or mixture—exactly the interior region of the Gamma bridge’s “power-law core + exponential cutoff.” (Clauset et al. 2009) So seeing a nearly straight log-log line is a clue worth probing the mechanism for, not an endpoint. Just how big a pitfall this “a clue is not a proof” is, the next section unfolds.

4.10 What counts as a “real” power law: fitting, testing, and pitfalls

If you need to judge whether real data (gene-expression tails, network degrees, …) is a power law, this section is required reading; if you only want the principles, skip to “Real-world cases.” The core reference is Clauset, Shalizi, Newman, Power-law distributions in empirical data (SIAM Review, 2009)—a methodological classic that threw cold water on the “everything is a power law” craze. (Clauset et al. 2009)

Two common practices are both wrong. A great many early “power-law” papers relied only on: (1) plotting the data on log-log and eyeballing a straight line; (2) doing least-squares linear regression on the log-log histogram and reporting a high \(R^2\). Neither holds up:

• log-log regression’s slope estimate is biased—binning artifacts, noisy tail bins, and the least-squares error assumptions simply do not hold in log-space;
• a high \(R^2\) cannot prove a power law—log-normal, exponential, and stretched-exponential all “look pretty straight” on log-log over a finite range. And the eye is even less of a test.

The correct four-step recipe (the current gold standard):

1. First fix the lower bound \(x_{\min}\): a power law usually holds only in the tail. For each candidate \(x_{\min}\), choose it by minimizing the Kolmogorov–Smirnov distance (data vs. fitted power law).
2. Estimate \(\lambda\) by maximum likelihood, not regression. In the continuous case this is the Hill estimator \(\hat\lambda = 1 + n\big/\sum_i \ln\frac{x_i}{x_{\min}}\), consistent and unbiased; regression is biased.
3. Goodness-of-fit test: compute the KS statistic, then obtain a \(p\)-value by bootstrap resampling. Small \(p\) → reject the power-law hypothesis. (Large \(p\) is only “not rejected,” not “proven to be a power law.”)
4. Model comparison: use likelihood-ratio tests to pit the power law against the log-normal, exponential, stretched exponential, and power law with exponential cutoff. The real question is not “is it a power law” but “is the power law better than these alternatives.”

The conclusion stings. Clauset et al. re-estimated 24 real data sets previously called power laws: only a few passed cleanly (word frequency being one of the best-supported); many were indistinguishable from log-normal, or fit better by a power law with exponential cutoff; a “pure power law” over the whole range is rare. The follow-up Broido & Clauset (2019) re-examined about a thousand networks, bluntly titled Scale-free networks are rare—genuine scale-free structure is uncommon. (Clauset et al. 2009; Broido and Clauset 2019)

Is it or isn’t it? The contrast is what’s interesting. Things that “look heavy-tailed” alike—some are power laws, some are not at all—and the criterion goes back to the process (the axis this chapter has been pushing):

Phenomenon	Often called	More defensible	Generating process
Word frequency (Zipf), earthquake energy	power law	power law (fairly well-supported)	criticality / preferential growth
The bulk of wealth, income	power law	log-normal (only the extreme upper tail nears Pareto)	multiplicative growth, no floor
Most network degree distributions	scale-free power law	often truncated or log-normal; pure scale-free rare	— (Broido–Clauset 2019)
Height, weight, measurement error	——	Gaussian	additive averaging → CLT
Radioactive lifetime, memoryless waiting	——	exponential	constant hazard rate

Four processes, four fates: adding many small quantities → Gaussian; memoryless waiting → exponential; multiplying many small quantities with no lower floor → log-normal; multiplicative + reset/floor (Kesten) or criticality/preferential attachment → a true power law. So the empirical question “is it a power law” is essentially “which process is behind it.” Heavy-tailed \(\ne\) power law: the log-normal is equally heavy-tailed yet comes from a completely different mechanism, and on finite data it is often hard to distinguish from a power law—this is exactly the core the Clauset test aims to pry apart, and why “seeing a line means probing the mechanism, not declaring victory.”

This does not say “power laws don’t exist,” but: evidence should use MLE + goodness-of-fit + model comparison, and honestly admit that alternative distributions often fit equally well or better. Note its two named prime rivals—log-normal and power law with exponential cutoff—are exactly the two already met in this chapter: the former from a multiplicative process + CLT (“Aside: wealth”), the latter the interior of the Gamma bridge. In other words, the empirical puzzle “is it a power law or one of those two” was already placed by this chapter’s maximum-entropy framework on the same family tree with all three. Practical tools: R’s poweRlaw package and Python’s powerlaw package both implement the recipe above.

The process table above is the “summary.” Next, “Real-world cases” makes two of its rows concrete—critical branching (producing a true power law \(T^{-3/2}\)) and the multiplicative + floor Kesten process (bulk log-normal, upper tail Pareto, often truncated)—to see how a concrete mechanism produces a power law, and incidentally how those two “rivals” emerge from the same batch of mechanisms.

4.11 Real-world case: how critical branching in cell lineages gives a power-law tail

This makes the first row (criticality / branching) of the previous section’s process table concrete: a power law can be derived not only from the maximum entropy of “knowing only the geometric scale,” but can also grow out of a concrete cell-division mechanism. This section’s conclusion in one sentence: an exactly critical (division = termination) branching process has a cumulative output following a \(T^{-3/2}\) power law; deviate slightly from criticality and the tail is exponentially truncated (landing in the Gamma bridge’s “power-law core + exponential tail”). This is a classical Galton–Watson branching-process result. (Harris 1963) To see how this conclusion is computed exactly from combinatorics, read on; if you only want to remember this sentence, skip to “Aside: wealth.”

Consider a deliberately simplified cell-lineage model: start from one founder cell; each cell independently divides into two daughters with probability \(q\), or terminates with probability \(1-q\) (interpretable as death or terminal differentiation); the two daughters then repeat the same rule.

What we study here is not the clone size at some instant, but the total number of cells ever produced until the lineage terminates, \(T\), counting the founder as one. If \(n\) divisions occur, the lineage has exactly \(n+1\) terminations, hence \(T=2n+1\). The number of binary lineage shapes satisfying this condition is the \(n\)-th Catalan number

\[ C_n=\frac{1}{n+1}\binom{2n}{n} \]

and each shape has the same probability. Therefore

\[ P(T=2n+1)=C_nq^n(1-q)^{n+1}. \]

The Catalan number is no mystery here: it counts all valid lineage trees in which “the living cells produced by division never run out before a termination event.” This closed form lets us see directly what happens at the critical point.

When \(q=1/2\), division and termination are exactly balanced—critical branching. Using \(C_n\sim4^n/(\sqrt{\pi}n^{3/2})\),

\[ P(T=2n+1)=C_n2^{-(2n+1)} \sim \sqrt{\frac{2}{\pi}}\,T^{-3/2}. \]

So this discrete mechanism produces an asymptotic tail with power exponent \(3/2\). It falls in the earlier \(1<\lambda\le2\) range: in an idealized, infinitely continuing critical model, \(E[T]\) does not exist. This does not mean a real piece of tissue produces infinitely many cells, but that a tiny fraction of extremely long lineages is enough to rob the theoretical mean of a stable scale.

## ── Cell lineage: critical branching gives a power law; deviation from criticality brings a cutoff ──

n_branch <- 0:1500
t_branch <- 2 * n_branch + 1
q_branch <- c(0.50, 0.45, 0.40)
cols_branch <- c("tomato", "steelblue", "gray40")

log_p_total <- function(n, q) {
  lchoose(2 * n, n) - log(n + 1) + n * log(q) + (n + 1) * log1p(-q)
}
p_branch <- sapply(q_branch, function(q) exp(log_p_total(n_branch, q)))

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.6, 1))

plot(t_branch, p_branch[, 1], type = "l", log = "xy", col = cols_branch[1], lwd = 3,
     xlim = c(1, max(t_branch)), ylim = c(1e-30, 1),
     main = "Total lineage output T",
     xlab = "Total cells ever produced (log scale)",
     ylab = "P(T = t), log scale",
     cex.lab = 1.18, cex.main = 1.05)
for (i in 2:length(q_branch)) {
  lines(t_branch, p_branch[, i], col = cols_branch[i], lwd = 2.7)
}
legend("topright", c("q = 0.50: critical", "q = 0.45: subcritical", "q = 0.40: subcritical"),
       col = cols_branch, lwd = c(3, 2.7, 2.7), bty = "n", cex = 0.95)

scaled_p <- sweep(p_branch, 1, t_branch^(3/2), "*")
plot(t_branch, scaled_p[, 1], type = "l", log = "x", col = cols_branch[1], lwd = 3,
     xlim = c(1, max(t_branch)), ylim = c(0, 0.9),
     main = "Remove the T^(-3/2) factor",
     xlab = "Total cells ever produced (log scale)",
     ylab = "T^(3/2) P(T = t)",
     cex.lab = 1.18, cex.main = 1.05)
for (i in 2:length(q_branch)) {
  lines(t_branch, scaled_p[, i], col = cols_branch[i], lwd = 2.7)
}
abline(h = sqrt(2 / pi), col = "tomato", lty = 2, lwd = 2)
text(25, 0.85, "critical limit = sqrt(2/pi)", col = "tomato", cex = 0.95)

In the left panel, the red line for \(q=1/2\) approaches a straight line on log-log; at \(q=0.45\) or \(0.40\) the tail bends downward. The right panel divides out the \(T^{-3/2}\) power-law factor: the critical curve approaches the constant \(\sqrt{2/\pi}\), while the two subcritical curves keep dropping to 0. This is not an eyeballed fit but the exact computation from the Catalan formula above.

More precisely, for \(q<1/2\),

\[ P(T=t)\asymp t^{-3/2}\exp\left[-\frac{t}{2}\log\frac{1}{4q(1-q)}\right]. \]

That is, “power-law core \(\times\) exponential cutoff.” The further \(q\) is from \(1/2\), the earlier the cutoff appears; if \(q>1/2\), there is positive probability of never terminating, and \(T\) is no longer a probability distribution taking only finite values. Real tissue is further limited by space, resources, developmental time windows, and regulatory feedback, so even if locally near critical, the observed tail is further truncated.

This case answers a different question than the earlier maximum-entropy derivation. The earlier continuous Pareto was: on \([x_{\min},\infty)\) with only \(E[\ln X]\) fixed, the distribution with the fewest extra assumptions; here it is: given a microscopic cell division/termination rule, what the cumulative lineage output looks like. The former gives an inference principle, the latter a cellular mechanism producing an approximate power law. Neither replaces the other.

4.11.1 Aside: why is wealth so often discussed as a power law?

This makes the second and third rows of the process table concrete (multiplicative growth, and multiplicative + floor). Wealth is a very useful contrast, correcting a common misstatement: “multiplicative growth necessarily produces a power law” is not true. With pure proportional growth

\[ W_{t+1}=A_tW_t, \]

then \(\ln W_t=\ln W_0+\sum_{s<t}\ln A_s\); when many small independent growth rates add, the central limit theorem points first to a log-normal, not a power law.

A classic class of mechanisms producing a Pareto tail adds, on top of proportional growth, a continuous inflow of income, a minimum guarantee, or re-entry:

\[ W_{t+1}=A_tW_t+B_t,\qquad B_t>0. \]

Under suitable conditions (e.g. \(E[\ln A]<0\) so the overall process has a stationary state, yet with some periods where \(A>1\)), the upper tail of this Kesten-type recursion satisfies

\[ P(W>w)\sim w^{-\alpha}, \qquad E[A^\alpha]=1, \]

where \(\alpha\) is the tail/survival exponent (corresponding to the \(\lambda-1\) of the earlier notation convention, i.e. density exponent \(\lambda=\alpha+1\)). Intuitively, \(B_t\) keeps wealth from falling to zero, while occasional runs of positive growth push a few trajectories very high; together they form a steady state that is “bulk with a scale, upper tail approximately scale-free.” (Kesten 1973; Gabaix 2009)

So what is common in wealth data is not a pure power law from bottom to top, but possibly a bulk that is approximately log-normal, nearing a Pareto tail only above a high threshold; taxation, bankruptcy, finite lifespan, and wealth ceilings all bend the tail down. Together with the cell-lineage example, it shows: a power law is not a universal label but a tail description jointly determined by the concrete mechanism, the observation scale, and the truncation conditions.

---

5 Normal distribution = maximum entropy + mean + variance constraint

5.1 Derivation

Constraints: \(E[X] = \mu\), \(\text{Var}(X) = \sigma^2\). Two constraint functions are needed, \(f_1(x) = x\) and \(f_2(x) = x^2\); plug into the general solution:

\[p(x) = \frac{1}{Z}\, e^{-\lambda_1 x - \lambda_2 x^2}\]

Completing the square: the exponent is a quadratic in \(x\); complete the square (fold the linear term into the square), absorbing the constant term into \(Z\):

\[-\lambda_1 x - \lambda_2 x^2 = -\lambda_2\Bigl(x + \tfrac{\lambda_1}{2\lambda_2}\Bigr)^2 + \underbrace{\tfrac{\lambda_1^2}{4\lambda_2}}_{\text{absorbed into }Z}\]

So \(p(x) \propto e^{-\lambda_2 (x - \mu)^2}\), a Gaussian centered at \(\mu = -\tfrac{\lambda_1}{2\lambda_2}\) with variance \(\tfrac{1}{2\lambda_2}\). Matching the constraints gives \(\lambda_2 = \tfrac{1}{2\sigma^2}\), \(\lambda_1 = -\tfrac{\mu}{\sigma^2}\), and normalization constant \(Z = \sqrt{2\pi\sigma^2}\):

\[p(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\, e^{-(x - \mu)^2 / (2\sigma^2)}\]

Differential entropy \(H = -\int p\ln p\,dx\) (the steps are isomorphic to 3.1, just replacing “mean” with “variance”):

\[ \begin{aligned} H &= -\int p\,\ln p\,dx \\[2pt] &= -\int p\Bigl(-\tfrac12\ln(2\pi\sigma^2) - \tfrac{(x-\mu)^2}{2\sigma^2}\Bigr)dx && \ln p = -\tfrac12\ln(2\pi\sigma^2) - \tfrac{(x-\mu)^2}{2\sigma^2} \\[2pt] &= \tfrac12\ln(2\pi\sigma^2)\!\int p\,dx + \tfrac{1}{2\sigma^2}\!\int (x-\mu)^2 p\,dx && \text{split up} \\[2pt] &= \tfrac12\ln(2\pi\sigma^2)\cdot 1 + \tfrac{1}{2\sigma^2}\cdot\sigma^2 && \textstyle\int p\,dx = 1,\ \int (x-\mu)^2 p\,dx = \sigma^2 \\[2pt] &= \tfrac12\ln(2\pi\sigma^2) + \tfrac12 = \tfrac12\ln(2\pi e\,\sigma^2) && \tfrac12 = \tfrac12\ln e \end{aligned} \]

That is, \(H[\text{Normal}(\mu,\sigma^2)] = \tfrac12\ln(2\pi e\,\sigma^2)\). Note it depends only on \(\sigma^2\) and not on \(\mu\)—translating the distribution does not change how spread out it is.

Watch next (optional). Why π is in the normal distribution (beyond integral tricks) (3Blue1Brown) explains where the \(\pi\) in the normalization constant above comes from; it fits best once the Gaussian form itself feels familiar.

5.2 Numerical check: among the candidates, the normal has the highest entropy

Compare entropies among all distributions with \(\mu = 0\), \(\sigma^2 = 1\):

x_full <- seq(-8, 8, length.out = 8000)
dx_f   <- diff(x_full)[1]

h_full <- function(p_vals) {
  p_vals <- pmax(p_vals, 1e-300)
  p_norm <- p_vals / (sum(p_vals) * dx_f)
  -sum(p_norm * log(p_norm)) * dx_f
}

# Various distributions with mean = 0, variance = 1
sigma_target <- 1

# Laplace: Var = 2b^2, so b = 1/sqrt(2)
b_lap <- 1 / sqrt(2)
# Logistic: Var = pi^2 * s^2 / 3, so s = sqrt(3)/pi
s_log <- sqrt(3) / pi
# Uniform: Var = (2a)^2/12, so a = sqrt(3)
a_unif <- sqrt(3)

dists_n <- list(
  "Normal(0,1)"    = dnorm(x_full, 0, sigma_target),
  "Laplace(0, b)"  = 0.5/b_lap * exp(-abs(x_full)/b_lap),
  "Logistic(0, s)" = dlogis(x_full, 0, s_log),
  "Uniform[-a, a]" = dunif(x_full, -a_unif, a_unif),
  "t(df=5) scaled" = dt(x_full * sqrt(5/3), df = 5) * sqrt(5/3)
)

ent_n    <- sapply(dists_n, h_full)
h_theory_n <- 0.5 * log(2 * pi * exp(1))

par(mfrow = c(1, 2), mar = c(4, 7, 3.5, 1))

plot(NA, xlim = c(-4, 4), ylim = c(0, 0.45),
     main = "Distributions with mean=0, var=1",
     xlab = "x", ylab = "Density")
dist_cols_n <- c("tomato", "steelblue", "forestgreen", "orange", "purple")
for (i in seq_along(dists_n)) {
  lines(x_full, dists_n[[i]], col = dist_cols_n[i], lwd = 2,
        lty = ifelse(i == 1, 1, 2))
}
legend("topright", names(dists_n), col = dist_cols_n, lwd = 2,
       lty = c(1, rep(2, 4)), bty = "n", cex = 0.95)

barplot(rev(ent_n), horiz = TRUE, col = rev(dist_cols_n),
        border = "white", las = 1, xlim = c(0, max(ent_n) * 1.1),
        main = "Entropy comparison (mean=0, var=1)",
        xlab = "Differential entropy H")
abline(v = h_theory_n, col = "tomato", lwd = 2, lty = 2)
text(h_theory_n * 1.03, 3.0, paste("Normal theoretical =", round(h_theory_n, 3)),
     col = "tomato", adj = 0.5, srt = 90, cex = 0.95)

The red bar for the normal is longest—with both mean and variance fixed, the normal distribution really does have the maximum entropy.

Laplace and Logistic, though close in shape, differ in tail behavior—either too heavy (Laplace) or too light (Uniform), both injecting extra assumptions.

5.3 CLT: another road to the same destination

The central limit theorem says: no matter the original distribution, the mean of many independent samples tends to the normal.

Watch next. But what is the Central Limit Theorem? (3Blue1Brown) animates why very different shapes converge under summation; watch it before the simulation below or use it afterward as a visual cross-check.

n_sim <- 5000

par(mfrow = c(2, 3), mar = c(4, 3, 3, 1))

for (gen_info in list(
  list(gen = function(n) runif(n, -1, 1),       name = "Uniform(-1, 1)"),
  list(gen = function(n) rexp(n) - 1,            name = "Exp(1) - 1  [right skew]"),
  list(gen = function(n) rbinom(n,1,0.1) - 0.1,  name = "Bernoulli(0.1)  [extreme skew]")
)) {
  for (n_avg in c(5, 50)) {
    mat <- matrix(gen_info$gen(n_sim * n_avg), nrow = n_sim, ncol = n_avg)
    z   <- scale(rowMeans(mat))
    hist(z, breaks = 55, freq = FALSE, col = "steelblue", border = "white",
         main = paste0(gen_info$name, "  (n=", n_avg, ")"),
         xlab = "Standardized mean", xlim = c(-4, 4), ylim = c(0, 0.47),
         cex.main = 0.95)
    curve(dnorm(x), add = TRUE, col = "tomato", lwd = 2)
  }
}

CLT and maximum entropy are two independent roads to the normal distribution:

• CLT says: the limit of addition is normal (a mathematical fact)
• maximum entropy says: at fixed mean and variance, the least-biased distribution is normal (a logical principle)

Their conclusions agree, but their logic is entirely independent. CLT needs no concept of entropy, and maximum entropy needs no concept of summation.

Below is a proof of the CLT via the Fourier transform. If you only want the conclusion, skip straight to the next section, “Physical examples.”

5.3.1 Details: proving the CLT via the Fourier transform

The idea. As said before: adding independent random variables convolves their densities. Adding \(n\) of them is the density convolved with itself \(n\) times—computing this directly is a nightmare. But the Fourier transform has a signature property: it turns convolution into multiplication. So “convolve \(n\) times” becomes “raise to the \(n\)-th power,” immediately tractable. This is the engine of the whole proof.

Watch first (optional). But what is a convolution? (3Blue1Brown) gives the geometric picture of convolution. It is not needed for the CLT conclusion above, but it makes the Fourier proof below easier to follow.

What to prove. \(X_1,\dots,X_n\) i.i.d. with mean \(\mu\) and variance \(\sigma^2\) (finite). The standardized sum \(Z_n = \dfrac{\sum_i X_i - n\mu}{\sigma\sqrt{n}}\) converges in distribution to \(N(0,1)\). First center and standardize \(Y_i = (X_i-\mu)/\sigma\) (\(E[Y_i]=0\), \(E[Y_i^2]=1\)), so \(Z_n = \tfrac{1}{\sqrt n}\sum_i Y_i\).

Tool: the Fourier transform of the density \(\hat f(t) = \int f(x)\,e^{itx}\,dx = E[e^{itX}]\) (in probability it is called the characteristic function, essentially just the Fourier transform of the density). We use two signature Fourier properties, plus a convergence bridge:

1. Convolution theorem: \(\widehat{f * g} = \hat f\cdot\hat g\) (convolution in the time domain = multiplication in the frequency domain). So the transform of an independent sum = the product of the individual transforms—this turns “sum” into “product.”
2. Gaussian self-transform: the Fourier transform of the \(N(0,1)\) density is \(e^{-t^2/2}\) (a Gaussian is still a Gaussian after transforming).
3. Lévy continuity theorem (the convergence bridge, treated as a black box): pointwise convergence of \(\hat f\) \(\Rightarrow\) convergence in distribution.

Computation. By the convolution theorem, \(\hat f_{Z_n}(t) = \Bigl[\hat f_Y\!\bigl(\tfrac{t}{\sqrt n}\bigr)\Bigr]^n\). Taylor-expand \(\hat f_Y\) at 0 (this step uses mean 0, variance 1):

\[ \begin{aligned} \hat f_Y(s) = E[e^{isY}] &= 1 + is\,E[Y] - \tfrac{s^2}{2}E[Y^2] + o(s^2) && \text{expand }e^{isY} \\[2pt] &= 1 - \tfrac{s^2}{2} + o(s^2) && E[Y]=0,\ E[Y^2]=1 \end{aligned} \]

Substituting \(s = t/\sqrt n\) gives \(\hat f_Y(\tfrac{t}{\sqrt n}) = 1 - \tfrac{t^2}{2n} + o(\tfrac1n)\), so

\[\hat f_{Z_n}(t) = \Bigl[1 - \frac{t^2}{2n} + o\!\bigl(\tfrac1n\bigr)\Bigr]^n \;\xrightarrow{\,n\to\infty\,}\; e^{-t^2/2}\]

(using \((1+\tfrac an+o(\tfrac1n))^n \to e^a\) with \(a=-t^2/2\).)

At this point we have only proved pointwise convergence of the transform \(\hat f_{Z_n}(t)\to e^{-t^2/2}\); one last step turns it into convergence in distribution. \(e^{-t^2/2}\) is exactly the transform of \(N(0,1)\) (property 2, responsible for “recognizing” which distribution it is); then by the Lévy continuity theorem (tool 3)—pointwise convergence of characteristic functions, with the limit continuous at \(t=0\) (\(e^0=1\) ✓), implies convergence in distribution—we get \(Z_n \xrightarrow{d} N(0,1)\). ∎ (Continuity at \(t=0\) is not decorative: it rules out the degenerate case of probability mass escaping to infinity.)

Why the Gaussian in particular. Look at the expansion \(\hat f_Y(s) = 1 - \tfrac{s^2}{2} + \cdots\): the \(\tfrac{1}{\sqrt n}\) scaling is tuned precisely so that the second-order term (the variance) exactly survives the \(n\)-th power and gathers into \(e^{-t^2/2}\), while third order and higher (skewness…) carry higher powers of \(1/\sqrt n\) and are crushed to zero. Only “mean 0, variance 1” survives in the limit, everything else erased—so no matter what the original distribution looks like (as long as the variance is finite), the sum converges to the same Gaussian. This dovetails with maximum entropy: fix mean and variance, ignore the rest \(\to\) normal. Two Easter eggs: (1) \(e^{-t^2/2}\) is itself a Gaussian, and the quadratic in the exponent is the same quadratic as in the normal’s \(\ln p\); (2) the Gaussian is a fixed point of convolution (two Gaussians convolve to a Gaussian), so the limit stably settles on the Gaussian and changes no further.

5.4 Physical examples: chemical-bond vibration and Brownian diffusion

Chemical-bond vibration. Approximate a chemical bond as a spring; the potential energy for a displacement \(x\) from equilibrium is quadratic: \(E = \frac{1}{2}\kappa x^2\). Immersed in a thermal bath at temperature \(T\), the Boltzmann factor acts on the quadratic energy:

\[p(x) \propto e^{-\kappa x^2 / 2k_BT} = \text{Normal}\!\left(0,\ \sigma^2 = k_BT/\kappa\right)\]

A molecular-scale instance of “quadratic constraint → normal”: atoms undergo Gaussian-distributed thermal fluctuations around the equilibrium bond length, with the distribution widening at higher temperature (\(\sigma \propto \sqrt{T}\)), directly measurable in vibrational spectra. The velocity components of gas molecules follow the same logic (the kinetic energy \(\frac{1}{2}mv_x^2\) is also quadratic), see the “Boltzmann distribution” chapter.

Why does this give a Gaussian rather than “the Boltzmann distribution”? Because the Boltzmann factor is the mechanism, and which named distribution you get is decided by the degree of \(E(x)\). All \(e^{-E/k_BT}\) says is that probability decays exponentially in energy; it says nothing about how the coordinate \(x\) is distributed—that depends on how energy depends on \(x\):

• potential linear in the coordinate (the isothermal atmosphere’s \(E = mgh\)) → \(e^{-mgh/k_BT}\) is an exponential distribution in \(h\);
• potential quadratic in the coordinate (here \(E = \frac12\kappa x^2\), or the kinetic energy \(\frac12 mv_x^2\)) → \(e^{-\kappa x^2/2k_BT}\) is a Gaussian in \(x\).

One factor, one substitution of \(E(x)\), a different named distribution. So the exponential and the Gaussian are not competitors of the Boltzmann distribution—they are two faces of it under different forms of the energy. The next chapter walks this line all the way through, Maxwell speed distribution included.

Brownian diffusion. Drop ink into water, and each pigment particle’s displacement is the accumulation of impulses from countless molecular collisions—precisely the physical incarnation of the CLT of the previous subsection (Einstein 1905). The position distribution at time \(t\) is a Gaussian with variance \(\propto t\), its profile widening as \(\sqrt{t}\).

## ── Physical examples: chemical-bond vibration + Brownian diffusion ──────────

par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))

# Left: harmonic-oscillator potential E = kappa x^2 / 2 in a thermal bath → position distribution is normal
kappa <- 1
x_b   <- seq(-3.5, 3.5, length.out = 400)
E_pot <- 0.5 * kappa * x_b^2
plot(x_b, E_pot / max(E_pot) * 0.9, type = "l", lwd = 2.5, col = "gray50",
     ylim = c(0, 1.0),
     main = "Bond vibration:  p(x) ~ exp(-kappa x^2 / 2kT)",
     xlab = "Displacement x from equilibrium", ylab = "(rescaled)",
     cex.lab = 1.35, cex.main = 1.1)
kT_b   <- c(0.3, 1.5)
cols_b <- c("steelblue", "tomato")
for (i in seq_along(kT_b)) {
  sd_i <- sqrt(kT_b[i] / kappa)
  lines(x_b, dnorm(x_b, 0, sd_i) / dnorm(0, 0, sqrt(kT_b[1]/kappa)) * 0.9,
        col = cols_b[i], lwd = 3)
}
legend("topright",
       c("potential E = kappa x^2/2",
         paste0("p(x) at kT = ", kT_b, "  (sd = ", round(sqrt(kT_b/kappa), 2), ")")),
       col = c("gray50", cols_b), lwd = c(2.5, 3, 3), bty = "n", cex = 1.0)

# Right: Brownian diffusion—the random-walk position distribution is a widening Gaussian
n_part  <- 20000
t_snap  <- c(25, 100, 400)
cols_d  <- c("steelblue", "purple", "tomato")
# Each particle = accumulation of many independent collision impulses (each step ~ N(0,1))
pos_at <- lapply(t_snap, function(t_i)
  rowSums(matrix(rnorm(n_part * t_i), nrow = n_part)))
plot(NA, xlim = c(-70, 70), ylim = c(0, 0.085),
     main = "Diffusion: Gaussian spreading, sd ~ sqrt(t)",
     xlab = "Position", ylab = "Density",
     cex.lab = 1.35, cex.main = 1.1)
for (i in seq_along(t_snap)) {
  d_i <- density(pos_at[[i]])
  lines(d_i, col = cols_d[i], lwd = 3)
  curve(dnorm(x, 0, sqrt(t_snap[i])), add = TRUE,
        col = cols_d[i], lwd = 1.5, lty = 2)
}
legend("topright",
       c(paste0("t = ", t_snap, "  (sd = ", sqrt(t_snap), ")"),
         "N(0, t) theory (dashed)"),
       col = c(cols_d, "gray40"), lwd = c(3, 3, 3, 1.5),
       lty = c(1, 1, 1, 2), bty = "n", cex = 1.0)

---

6 Boltzmann distribution = maximum entropy + energy constraint

6.1 The Boltzmann principle is the maximum entropy principle

A thermodynamic system at temperature \(T\):

\[\boxed{P(\text{state}) \propto e^{-E / k_B T}}\]

This formula is often memorized as a “physical law.” But its derivation is just maximum entropy:

• constraint: total energy conserved, \(\langle E \rangle\) fixed
• maximum-entropy solution: \(p \propto e^{-\lambda E}\)

First, how this chapter sits relative to the earlier ones. In the opening table, the first three rows constrain an explicit function of \(x\) (\(x\), \(\ln x\), \(x^2\)), whereas this row reads \(f(x) = E(x)\)—a master statement, not a fourth named distribution standing beside the others. Once \(E(x)\) is specified, what falls out is one of the distributions already met:

How energy depends on the variable	\(e^{-E/k_BT}\) becomes	Instance
\(E\) linear in the variable: \(E = mgh\), evenly spaced levels	an exponential distribution in that variable	isothermal atmosphere, level occupancy
\(E\) quadratic: \(E = \frac12 mv_x^2\), \(\frac12\kappa x^2\)	a Gaussian in that variable	velocity components, bond vibration
quadratic energy \(+\) non-uniform density of states \(g(v)\propto v^2\)	the Maxwell speed distribution	molecular speeds

So “exponential, Gaussian, Boltzmann” are not three competing distributions: the Boltzmann factor is the shared mechanism, and the exponential and the Gaussian are its two faces under linear and quadratic energy, with the Maxwell speed distribution arising once the density of states is folded in. This chapter walks through them in that order.

Following that thread, one thing deserves to be said plainly: the Boltzmann form by itself rules out no distribution at all. Given any strictly positive density \(p(x)\), simply define

\[E(x) \;\equiv\; -k_BT\ln p(x)\]

and immediately \(p(x)\propto e^{-E(x)/k_BT}\). In other words, “this is a Boltzmann distribution,” read as a claim about the variable \(x\), is a tautology—every positive density can be written that way.

That makes it a different kind of object from a family such as the Gamma:

	Type	How it “contains” other distributions	Strength as a claim
Gamma	a genuine two-parameter family \(x^{k-1}e^{-\lambda x}\)	turn the parameters: \(k=1\) gives the exponential, \(k=\nu/2\) the chi-squared (note it does not contain the normal; \(k\to\infty\) is only asymptotic)	a strong claim, falsifiable by data
Boltzmann	a recipe with a free function \(E(x)\)	swap the function: \(E\) linear gives the exponential, quadratic gives the normal (both exactly)	with \(E\) unspecified, rules out nothing

The whole force of the Boltzmann factor therefore lies not in the exponential form but in physics telling you, independently, what \(E(x)\) is—gravity supplies \(mgh\), a spring supplies \(\frac12\kappa x^2\), kinetic energy supplies \(\frac12 mv^2\). \(E\) is an external input, not something fitted to the data; once \(E\) is pinned down, \(e^{-E/k_BT}\) becomes a distribution with temperature as its only free parameter, and it can be falsified.

The real umbrella, then, is not the Boltzmann distribution but the exponential family of the closing chapter, \(p\propto m(x)e^{-\sum_i\lambda_i f_i(x)}\): the exponential (\(f=x\)), power law (\(f=\ln x\)), normal (\(f=x,x^2\)) and Gamma (\(f=x,\ln x\)) are siblings under that umbrella—each a different choice of \(f\)—while Boltzmann is the branch that leaves \(f\) blank and lets physics fill it in.

\(\lambda\) is the Lagrange multiplier. Physicists gave it a name: \(\beta = 1/(k_BT)\).

That is: temperature is not an independent physical quantity but the reciprocal of the Lagrange multiplier in the maximum-entropy derivation.

• high temperature (\(\beta\) small) → loose constraint → more spread-out distribution → closer to equiprobable
• low temperature (\(\beta\) large) → tight constraint → distribution concentrated in low-energy states
• \(T \to \infty\) → \(\beta \to 0\) → \(p \to \text{Uniform}\) (all states equiprobable)

par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))

E_seq   <- seq(0, 10, length.out = 400)
kT_vals <- c(0.5, 1, 2, 4)
cols    <- c("steelblue", "tomato", "forestgreen", "purple")

# Left: energy distributions at different temperatures
plot(NA, xlim = c(0, 10), ylim = c(0, 2.1),
     main = "Boltzmann: P(E) = beta * exp(-beta * E)\nbeta = Lagrange multiplier = 1/kT",
     xlab = "Energy E", ylab = "P(E)")
for (i in seq_along(kT_vals)) {
  lines(E_seq, dexp(E_seq, rate = 1/kT_vals[i]), col = cols[i], lwd = 2.5)
}
legend("topright", paste("kT =", kT_vals, "  (beta =", 1/kT_vals, ")"),
       col = cols, lwd = 2.5, bty = "n", cex = 0.95)

# Right: temperature vs entropy
# H[Exp(beta)] = 1 + ln(1/beta) = 1 + ln(kT)
kT_range <- seq(0.2, 5, length.out = 200)
H_boltz  <- 1 + log(kT_range)

plot(kT_range, H_boltz, type = "l", col = "tomato", lwd = 2.5,
     main = "Temperature controls entropy\nH = 1 + ln(kT)",
     xlab = "kT  (temperature)", ylab = "Entropy H")
abline(h = 0, col = "gray70", lty = 3)
text(3.5, 0.5, "Higher T = higher entropy\n= more \"spread out\"",
     col = "gray40", cex = 0.95)

6.2 What is temperature? \(1/T = \partial S / \partial E\)

“\(\lambda\) is just \(1/k_BT\)” is not merely a change of notation—it is the strict identity between the thermodynamic definition of temperature and the maximum-entropy multiplier. Compute the entropy directly for the maximum-entropy solution \(p(E) = \frac{1}{Z}e^{-\lambda E}\):

\[H = -\int p \ln p = -\int p\,\bigl(-\lambda E - \ln Z\bigr) = \lambda \langle E\rangle + \ln Z\]

Differentiate with respect to \(\langle E\rangle\) (note \(Z\) depends on \(\langle E\rangle\) through \(\lambda\), using \(\frac{d\ln Z}{d\lambda} = -\langle E\rangle\)):

\[\frac{dH}{d\langle E\rangle} = \lambda + \underbrace{\left(\langle E\rangle + \frac{d\ln Z}{d\lambda}\right)}_{=\,0}\frac{d\lambda}{d\langle E\rangle} = \lambda\]

And the thermodynamic definition of temperature (Clausius, \(dS = \delta Q / T\)) is exactly \(\dfrac{1}{T} = \dfrac{\partial S}{\partial E}\). Matching the two (with \(S = k_B H\)):

\[\lambda = \frac{1}{k_B T}\]

Verify with this section’s exponential distribution: \(H = 1 + \ln\langle E\rangle\), \(\frac{dH}{d\langle E\rangle} = \frac{1}{\langle E\rangle} = \frac{1}{k_BT}\). ✓

Physical intuition: temperature measures “how much entropy rises per unit of energy injected”—a cold system rises a lot (\(1/T\) large), a hot system rises little. When two systems are in contact, energy flows from hot to cold because the same energy buys more entropy on the cold side, so the total entropy increases—the direction of heat conduction is also derived from maximum entropy, needing no extra assumption.

6.3 Linear energy levels → exponential distribution

When the density of energy states is uniform (e.g. the equally spaced levels of a quantum harmonic oscillator, or the kinetic energy of a 2D gas), after normalization:

\[P(E) = \beta\, e^{-\beta E} = \text{Exponential}(\beta = 1/k_BT)\]

This is the exponential distribution. The Boltzmann distribution and the exponential distribution are not “similar”—they are the same thing.

kT <- 1

par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))

E_samples <- -kT * log(1 - runif(N))

hist(E_samples, breaks = 80, freq = FALSE, col = "steelblue", border = "white",
     main = "Sampling Boltzmann via  -kT * log(1 - U)",
     xlab = "Energy E")
curve(dexp(x, rate = 1/kT), add = TRUE, col = "tomato", lwd = 2)
legend("topright", "Exp(1/kT) theory", col = "tomato", lwd = 2, bty = "n")

# Arrhenius: reaction rate is a Boltzmann tail
Ea_vals <- c(0.5, 1, 2, 3)
kT_range2 <- seq(0.3, 5, length.out = 200)
plot(NA, xlim = c(0, 5), ylim = c(0, 1.05),
     main = "Arrhenius:  k ~ exp(-Ea/kT)\nReaction rate = Boltzmann tail",
     xlab = "kT  (temperature)", ylab = "exp(-Ea / kT)")
arr_cols <- c("steelblue","tomato","forestgreen","purple")
for (i in seq_along(Ea_vals)) {
  lines(kT_range2, exp(-Ea_vals[i] / kT_range2), col = arr_cols[i], lwd = 2)
}
legend("bottomright", paste("Ea =", Ea_vals), col = arr_cols, lwd = 2, bty = "n")

The Arrhenius formula for chemical reactions, \(k = A\, e^{-E_a/(k_BT)}\), is just the tail probability of the Boltzmann distribution—the fraction of molecules able to cross the energy barrier \(E_a\).

6.4 Quadratic kinetic energy → normal velocity distribution

The kinetic energy of a 3D gas molecule is a quadratic function of velocity: \(E_x = \frac{1}{2}mv_x^2\).

Applying the Boltzmann factor to the kinetic energy:

\[P(v_x) \propto e^{-mv_x^2 / (2k_BT)} = e^{-v_x^2 / (2\sigma^2)}, \quad \sigma = \sqrt{k_BT/m}\]

This is exactly the normal distribution. Compare the maximum-entropy view: \(f(x) = x^2\), Lagrange multiplier \(\lambda_2 = m/(2k_BT) = 1/(2\sigma^2)\)—the physics’ \(m/(2k_BT)\) is exactly the \(\lambda_2\) in the maximum-entropy derivation.

sigma <- 1
vx <- rnorm(N, 0, sigma)
vy <- rnorm(N, 0, sigma)
vz <- rnorm(N, 0, sigma)

par(mfrow = c(1, 3), mar = c(4, 4, 3.5, 1))
for (comp in list(list(v = vx, name = "vx"),
                  list(v = vy, name = "vy"),
                  list(v = vz, name = "vz"))) {
  hist(comp$v, breaks = 80, freq = FALSE, col = "steelblue", border = "white",
       main = paste("Velocity", comp$name, "~ N(0, kT/m)"),
       xlab = comp$name, xlim = c(-4, 4))
  curve(dnorm(x, 0, sigma), add = TRUE, col = "tomato", lwd = 2)
}

6.5 Maxwell-Boltzmann speed distribution

The speed \(v = \sqrt{v_x^2 + v_y^2 + v_z^2}\), composed of three independent normal components:

\[f(v) = \sqrt{\frac{2}{\pi}} \frac{v^2}{\sigma^3} e^{-v^2/(2\sigma^2)}\]

speed  <- sqrt(vx^2 + vy^2 + vz^2)
mb_pdf <- function(v, sigma) sqrt(2/pi) * v^2 / sigma^3 * exp(-v^2 / (2*sigma^2))

v_peak <- sqrt(2) * sigma
v_mean <- 2/sqrt(pi) * sigma
v_rms  <- sqrt(3) * sigma

par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))

v_seq <- seq(0, 5.5, length.out = 400)
hist(speed, breaks = 80, freq = FALSE, col = "steelblue", border = "white",
     main = "Maxwell-Boltzmann speed distribution",
     xlab = "Speed v", xlim = c(0, 5.5))
lines(v_seq, mb_pdf(v_seq, sigma), col = "tomato", lwd = 2.5)
abline(v = v_peak, col = "orange",      lty = 2, lwd = 1.8)
abline(v = v_mean, col = "forestgreen", lty = 2, lwd = 1.8)
abline(v = v_rms,  col = "purple",      lty = 2, lwd = 1.8)
legend("topright",
       c("MB theory",
         paste0("Mode  = sqrt(2)*sigma = ", round(v_peak, 2)),
         paste0("Mean  = 2*sigma/sqrt(pi) = ", round(v_mean, 2)),
         paste0("RMS   = sqrt(3)*sigma = ", round(v_rms, 2))),
       col = c("tomato","orange","forestgreen","purple"),
       lwd = c(2.5,1.8,1.8,1.8), lty = c(1,2,2,2), bty = "n", cex = 0.95)

plot(NA, xlim = c(0, 8), ylim = c(0, 0.85),
     main = "Effect of temperature on MB distribution",
     xlab = "Speed v", ylab = "Density")
kTm_vals <- c(0.5, 1, 2, 4)
for (i in seq_along(kTm_vals)) {
  lines(v_seq, mb_pdf(v_seq, sqrt(kTm_vals[i])), col = cols[i], lwd = 2.5)
}
legend("topright", paste("kT/m =", kTm_vals), col = cols, lwd = 2.5, bty = "n")

---

6.6 The partition function: not merely a normalizing constant

So far \(Z\) has been treated as “the denominator that makes the probabilities sum to one.” It is in fact the hardest-working tool in this chapter.

\[Z(\beta) \;=\; \sum_i e^{-\beta E_i} \qquad\text{(an integral in the continuous case)}\]

Differentiate \(\ln Z\) with respect to \(\beta\) and the mean energy falls out by itself:

\[-\frac{\partial \ln Z}{\partial \beta} = \frac{1}{Z}\sum_i E_i\,e^{-\beta E_i} = \bar E\]

Differentiate once more and out comes the fluctuation:

\[\frac{\partial^2 \ln Z}{\partial \beta^2} \;=\; \overline{E^2} - \bar E^2 \;=\; \mathrm{Var}(E)\]

One function, \(\ln Z\): one derivative gives the mean, two give the variance. That is no accident—it is the moment-generating machinery: \(\ln Z(\beta)\) is the cumulant generating function of the energy.

Try it on the smallest possible system: a two-level system (energies \(0\) and \(\Delta\)—say a molecule with only a ground and an excited state). Here \(Z = 1 + e^{-\beta\Delta}\), so

\[p_1 = \frac{e^{-\beta\Delta}}{1+e^{-\beta\Delta}}, \qquad \bar E = \Delta\,p_1\]

D <- 1                                    # gap Δ, with k_B = 1
Zf <- function(b) 1 + exp(-b*D)
p1 <- function(Tt) exp(-D/Tt)/Zf(1/Tt)
Em <- function(Tt) D*exp(-D/Tt)/Zf(1/Tt)
Va <- function(Tt) D^2*exp(-D/Tt)/Zf(1/Tt)^2
Ts <- seq(0.05, 3, length.out = 1200)
h  <- 1e-6
CV <- sapply(Ts, function(t) (Em(t+h) - Em(t-h))/(2*h))     # C_V = d<E>/dT, numerically

par(mfrow = c(1, 2), mar = c(4.6, 4.8, 3.6, 1))

plot(Ts, 1 - p1(Ts), type = "l", lwd = 3, col = "#2E86AB", ylim = c(0, 1.05),
     main = expression("Two-level system:  "*Z == 1 + e^{-Delta/kT}),
     xlab = expression("temperature  "*kT/Delta), ylab = "level occupancy",
     cex.lab = 1.2, cex.main = 1.1)
lines(Ts, p1(Ts), lwd = 3, col = "tomato")
abline(h = 0.5, col = "gray60", lty = 3)
text(2.35, 0.455, "both approach 1/2", col = "gray35", cex = 0.95)
text(0.95, 0.05, "frozen: no thermal budget", col = "gray35", cex = 0.95)
text(1.95, 0.30, expression(bar(E) == Delta%.%p[1]), col = "gray25", cex = 1.15)
legend(0.05, 0.92, legend = expression(p[0]~"(ground)", p[1]~"(excited)"),
       col = c("#2E86AB", "tomato"), lwd = 3, bty = "n", cex = 1.05)

plot(Ts, CV, type = "l", lwd = 4, col = "forestgreen", ylim = c(0, max(CV)*1.32),
     main = expression("Fluctuation = response:  "*Var(E) == k[B]*T^2*C[V]),
     xlab = expression("temperature  "*kT/Delta), ylab = expression(C[V]*"   or   "*Var(E)/T^2),
     cex.lab = 1.2, cex.main = 1.1)
lines(Ts, Va(Ts)/Ts^2, lwd = 2, col = "purple", lty = 2)
pk <- Ts[which.max(CV)]
points(pk, max(CV), pch = 16, col = "forestgreen", cex = 1.6)
arrows(pk + 0.60, max(CV)*1.13, pk + 0.06, max(CV)*1.02, length = 0.08, col = "gray45")
text(pk + 0.62, max(CV)*1.16, sprintf("Schottky peak\nkT = %.3f Delta", pk),
     col = "gray30", cex = 0.95, adj = 0)
legend("right", legend = expression(C[V] == d*bar(E)/dT, Var(E)/T^2),
       col = c("forestgreen", "purple"), lwd = c(4, 2), lty = c(1, 2), bty = "n", cex = 1.05)

Left panel: at low temperature the system is frozen into the ground state (\(p_1\to0\)—no thermal budget for climbing the gap); at high temperature the two levels are equally occupied (\(p_1\to1/2\), since \(k_BT \gg \Delta\) makes the gap irrelevant). Note that the high-temperature limit is not “everything piles into the upper level” but equiprobability—precisely the Boltzmann factor flattening out as \(T\to\infty\).

The right panel is the point of this section. The heat capacity \(C_V = \mathrm{d}\bar E/\mathrm{d}T\) is a macroscopic response you can measure: how much energy it takes to raise the temperature by one degree. \(\mathrm{Var}(E)\) is a microscopic fluctuation. The two curves lie exactly on top of each other because

\[\boxed{\;\mathrm{Var}(E) \;=\; k_B T^2\, C_V\;}\]

—fluctuation and response are two readings of the same second derivative, \(\partial^2\ln Z/\partial\beta^2\). A system that is easier to heat is one whose energy was already jittering more. The peak sits at \(k_BT \approx 0.417\,\Delta\) (the Schottky peak): too cold and nothing can climb, too hot and the levels are already saturated; only near the scale of the gap does a little added heat move the occupancies the most.

This “differentiate \(\ln Z\) to get moments” machine is not the property of thermodynamics. Remember its three parts—the \(\beta\) in the exponent, the \(E\) being averaged, the \(\ln Z\) that normalizes—because the closing chapter will show that statistics has its own names for all three (natural parameter, sufficient statistic, cumulant generating function), and that the Boltzmann distribution is simply the special case physics discovered first.

6.7 Degeneracy and density of states: the Boltzmann factor is only half of it

There is a gap here that is easy to miss. The Boltzmann factor says that each individual state has probability \(\propto e^{-\beta E}\). But if what you care about is not states but energy, you must also count: how many states share that energy.

\[p(E) \;\propto\; \underbrace{g(E)}_{\text{how many slots}}\;\cdot\;\underbrace{e^{-E/k_BT}}_{\text{weight of each}}\]

\(g(E)\) is the degeneracy (discrete levels) or the density of states (continuous energy). The two factors pull in opposite directions—more slots as \(E\) rises, less weight—so the peak of the probability is not at the lowest energy:

Tt <- 1; b <- 1/Tt; s <- 2                     # g(E) ∝ E^s, with k_B = 1
E  <- seq(0, 12, length.out = 1500)
g  <- E^s; bz <- exp(-b*E); pr <- g*bz
par(mfrow = c(1, 3), mar = c(4.5, 4.7, 3.6, 1))

plot(E, g/max(g), type = "l", lwd = 4, col = "#2E86AB",
     main = expression("How many states:  "*g(E)*" ~ "*E^2),
     xlab = "energy E", ylab = "g(E) (rescaled)", cex.lab = 1.25, cex.main = 1.1)
text(4.6, 0.72, "more room\nat high E", col = "gray30", cex = 1.05)

plot(E, bz, type = "l", lwd = 4, col = "tomato",
     main = expression("How likely each:  "*e^{-E/kT}),
     xlab = "energy E", ylab = "Boltzmann factor", cex.lab = 1.25, cex.main = 1.1)
text(6.2, 0.55, "penalty grows\nwith E", col = "gray30", cex = 1.05)

plot(E, pr/max(pr), type = "l", lwd = 4, col = "forestgreen",
     main = expression("Product:  "*p(E) %prop% g(E)*e^{-E/kT}),
     xlab = "energy E", ylab = "p(E) (rescaled)", cex.lab = 1.25, cex.main = 1.1)
polygon(c(E, rev(E)), c(pr/max(pr), rep(0, length(E))),
        col = adjustcolor("forestgreen", 0.18), border = NA)
abline(v = s*Tt, col = "gray40", lty = 2)
points(s*Tt, 1, pch = 16, col = "forestgreen", cex = 1.7)
text(s*Tt + 0.4, 0.90, expression("peak at  "*E^"*" == s*kT), col = "gray25", cex = 1.1, adj = 0)
text(4.6, 0.35, "NOT at E = 0:\ncompetition between\nroom and penalty",
     col = "gray30", cex = 1.02, adj = 0)

Substituting \(g(E)\propto E^s\) and maximizing gives the peak in closed form:

\[\frac{\mathrm{d}}{\mathrm{d}E}\Bigl[s\ln E - \beta E\Bigr] = 0 \;\Longrightarrow\; E^\ast = \frac{s}{\beta} = s\,k_BT\]

With \(s=2\) and \(k_BT=1\) the peak lands exactly at \(E^\ast=2\), as plotted. So the “most likely energy” rises linearly with temperature and is never zero—even though every individual state is more likely the lower it sits. This is where the \(v^2\) factor in the Maxwell speed distribution earlier in this chapter comes from: in three-dimensional velocity space the shell of radius \(v\) has area \(\propto v^2\), and that is \(g(v)\). Velocity components are pure Gaussians (uniform \(g\)); the speed has a peak.

Note carefully where this \(g\) sits: it comes neither from a constraint nor from the Boltzmann factor, but from what you choose as the unit of counting. The closing chapter gives it a name in statistics (the base measure) and shows that it has been quietly present in every derivation so far.

6.8 Equipartition: \(\frac12 k_BT\) per quadratic degree of freedom

Chapter 5 established that quadratic energy gives a Gaussian. Turn the question around—how much energy does such a degree of freedom hold on average?

For a one-dimensional spring \(E = \frac12\kappa x^2\), the position is Gaussian with \(\sigma^2 = k_BT/\kappa\), so

\[\bar E = \tfrac12\kappa\,\overline{x^2} = \tfrac12\kappa\sigma^2 = \tfrac12 k_BT\]

\(\kappa\) cancels. Whether the spring is stiff or soft, the molecule heavy or light, the energy held by that degree of freedom depends on temperature alone.

kap <- 2.7                                  # any stiffness, with k_B = 1
chk <- sapply(c(0.5, 1, 3), function(Tt){
  x <- rnorm(2e6, 0, sqrt(Tt/kap))          # p(x) ∝ exp(-kappa x^2 / 2T)
  c(T = Tt, simulated = mean(0.5*kap*x^2), theory = 0.5*Tt)
})
round(t(chk), 5)

##        T simulated theory
## [1,] 0.5   0.24974   0.25
## [2,] 1.0   0.50003   0.50
## [3,] 3.0   1.49939   1.50

Generalized, this is the equipartition theorem: every quadratic term in the energy (each direction of kinetic energy, each harmonic coordinate) receives \(\frac12 k_BT\) on average, so \(f\) quadratic degrees of freedom give \(\bar E = \frac f2 k_BT\) and \(C_V = \frac f2 k_B\). A monatomic gas has three translational directions → \(\frac32 k_BT\), the familiar ideal-gas result.

It also explains something in passing: why temperature can be defined as “energy per degree of freedom.” \(k_BT\) is not an arbitrary conversion factor—it is the allowance the heat bath pays out to each quadratic degree of freedom.

6.9 Free energy: the other face of maximum entropy

One last step up in abstraction—it is this chapter’s bridge to the closing one.

Combine \(\bar E = -\partial_\beta\ln Z\) with the entropy (\(S = k_B H\), \(H = -\sum p\ln p\), substituting \(p_i = e^{-\beta E_i}/Z\)):

\[S = k_B\bigl(\beta\bar E + \ln Z\bigr) \qquad\Longrightarrow\qquad \boxed{\;-k_BT\ln Z \;=\; \bar E - TS \;\equiv\; F\;}\]

\(F\) is the free energy. The same fact now has two statements:

• maximum-entropy view: fix the mean energy \(\bar E\), maximize the entropy \(S\);
• free-energy view: fix the temperature \(T\), minimize \(F = \bar E - TS\).

The two are a Legendre transform pair (\(\bar E\) and \(\beta\) are conjugate variables, as price and quantity are). That is what lies behind the familiar physics slogan “systems minimize their free energy”—it is not a second law standing beside the maximum-entropy principle but the same principle rewritten for the case where temperature, rather than energy, is what is held fixed.

The minus sign in \(F\) acquires a meaning too: \(-TS\) is the subsidy paid to entropy. The hotter it gets, the more diversity is worth; as \(T\to0\) the entropy term loses force and only the lowest energy counts. This is the dual face of the statement in chapter 1 that temperature measures how much entropy grows per unit of energy injected.

The closing chapter names \(\ln Z\) in its statistical guise and shows how this Legendre duality becomes the convex conjugate of the exponential family—one structure, named once by each field.

7 The unified picture: the exponential family and \(p \propto e^{-\lambda f(x)}\)

7.1 One equation, four distributions

Back to the maximum-entropy general solution \(p \propto e^{-\lambda f(x)}\):

\[p(x) \propto e^{-\lambda_1 f_1(x) - \lambda_2 f_2(x) - \cdots}\]

In statistics this is called the exponential family. Nature reproduces these distributions repeatedly not because they are “pretty” but because they are the most impartial choice under limited information.

par(mar = c(1, 1, 2, 1))
plot(NA, xlim = c(0, 12), ylim = c(0, 8.5), axes = FALSE, bty = "n",
     xlab = "", ylab = "",
     main = "Unified view: MaxEnt + constraint => distribution")

draw_box <- function(x, y, label, col, w = 3.5, h = 0.7, cex = 0.95) {
  rect(x-w/2, y-h/2, x+w/2, y+h/2, col = col, border = "white", lwd = 2)
  text(x, y, label, col = "white", font = 2, cex = cex)
}
draw_arr <- function(x0, y0, x1, y1, label = "") {
  arrows(x0, y0, x1, y1, length = 0.12, lwd = 1.8, col = "gray40")
  mx <- (x0+x1)/2; my <- (y0+y1)/2
  if (nchar(label) > 0) text(mx + 0.15, my, label, cex = 0.95, col = "gray30", font = 3)
}

# Top: MaxEnt principle
draw_box(6, 8, "Maximize  H = -int p ln p  dx\nsubject to constraints",
         "gray25", w = 7.5, h = 0.9)

# General solution
draw_box(6, 6.3, "General solution:  p(x) ~ exp( -lambda * f(x) )",
         "gray45", w = 8, h = 0.7)

draw_arr(6, 7.55, 6, 6.7, "Lagrange multipliers")

# Four branches
draw_box(1.5, 4.3, "f(x) = x\nlambda = 1/mu",              "steelblue",  w = 2.55, cex = 0.95)
draw_box(4.5, 4.3, "f(x) = log x\nlambda > 1",             "purple",     w = 2.55, cex = 0.95)
draw_box(7.5, 4.3, "f(x) = x^2\nlambda = 1/(2*sigma^2)",    "tomato",     w = 2.55, cex = 0.95)
draw_box(10.5, 4.3, "f(x) = Energy\nlambda = beta = 1/kT",  "forestgreen", w = 2.55, cex = 0.95)

draw_arr(3.2, 5.95, 1.5,  4.65)
draw_arr(5.0, 5.95, 4.5,  4.65)
draw_arr(7.0, 5.95, 7.5,  4.65)
draw_arr(8.8, 5.95, 10.5, 4.65)

# Results
draw_box(1.5,  2.5, "Exponential(1/mu)",     "steelblue",  w = 2.55, cex = 0.95)
draw_box(4.5,  2.5, "Power law / Pareto",    "purple",     w = 2.55, cex = 0.95)
draw_box(7.5,  2.5, "Normal(mu, sigma^2)",   "tomato",     w = 2.55, cex = 0.95)
draw_box(10.5, 2.5, "Boltzmann = Exp(beta)", "forestgreen", w = 2.55, cex = 0.95)

draw_arr(1.5,  3.95, 1.5,  2.85)
draw_arr(4.5,  3.95, 4.5,  2.85)
draw_arr(7.5,  3.95, 7.5,  2.85)
draw_arr(10.5, 3.95, 10.5, 2.85)

# Physical examples
text(1.5,  1.42, "Waiting times\nDecay processes\nPoisson process",
     cex = 0.95, col = "gray40")
text(4.5,  1.42, "City sizes\nWealth tails\nNetwork degrees",
     cex = 0.95, col = "gray40")
text(7.5,  1.42, "Measurement error\nCLT (sum of many)\nDiffusion",
     cex = 0.95, col = "gray40")
text(10.5, 1.42, "Energy at temperature T\nArrhenius rates\nVelocity components",
     cex = 0.95, col = "gray40")

# Connection between Boltzmann → Exp and Normal
draw_arr(10.5, 2.15, 7.9, 2.15)
draw_arr(10.5, 2.05, 1.9, 2.05)
text(6.0, 0.5, "linear E => Exp;  quadratic E => Normal;  log-scale constraint => Power law",
     cex = 0.95, col = "gray50", font = 3)

7.2 First, the ruler: what a measure is

The next section needs the notion of a reference measure. The phrase sounds like a prerequisite from functional analysis; the object itself is far humbler—but leave it unexplained and the next section turns into an incantation. So this section makes it plain.

A measure is just “how you count how much there is.” Hand someone an interval \([a,b]\) and ask how big it is, and the natural answer is its length, \(b-a\). The person who made that idea of “assigning a size to a set” rigorous was Henri Lebesgue (1902), so the ruler that “counts by ordinary length” is called the Lebesgue measure—the unassuming \(\mathrm{d}x\) in every integral.

The point is this: there is more than one ruler. You can count by length, or by orders of magnitude (one notch per factor of ten), or one notch per integer. Change the ruler and the meaning of words like “uniform,” “spread out,” and “knowing nothing” changes with it.

7.2.1 Why entropy has to name its ruler

Discrete entropy never has to worry about this: it counts states, it is dimensionless, and no change of units touches it. The continuous case is different—differential entropy is not invariant under a change of coordinates. Watch:

Hnum <- function(dens, lo, hi)                     # differential entropy, numerically
  -integrate(function(x){p <- dens(x); ifelse(p > 0, p*log(p), 0)}, lo, hi)$value

# One and the same "uniformly random length between 0 and 1 metre", relabelled
c(in_metres      = Hnum(function(x) dunif(x, 0, 1),    0, 1),
  in_centimetres = Hnum(function(x) dunif(x, 0, 100),  0, 100),
  in_millimetres = Hnum(function(x) dunif(x, 0, 1000), 0, 1000),
  ln100 = log(100), ln1000 = log(1000)) |> round(4)

##      in_metres in_centimetres in_millimetres          ln100         ln1000 
##         0.0000         4.6052         6.9078         4.6052         6.9078

The same random fact, recorded in different units, has a different entropy—\(0\) in metres, \(\ln 100 = 4.605\) in centimetres. The reason is no mystery: a probability density carries units (probability per metre), so changing units rescales \(p\) and shifts \(\ln p\) bodily.

So the absolute value of a differential entropy means nothing by itself; it means something only once a ruler is fixed. That also explains an oddity met earlier: the exponential distribution has \(H<0\) when \(\mu < 1/e\). The minus sign does not mean “negative information”—it means this distribution is more concentrated than “one unit of length wide.” Differential entropy measures spread relative to the ruler.

7.2.2 Relative entropy is the quantity that does not depend on the ruler

The fix is to write the ruler into the formula. Relative entropy (KL divergence) does exactly that:

\[D(p\,\|\,m) \;=\; \int p(x)\,\ln\frac{p(x)}{m(x)}\,\mathrm{d}x\]

Here \(m\) is the ruler (the reference measure). Under a change of units \(p\) and \(m\) rescale together, the ratio \(p/m\) is untouched, and so is \(D\):

# p = Beta(2,5); reference m = uniform on the same interval. Then redo it in centimetres
D_m  <- integrate(function(x) dbeta(x,2,5)*log(dbeta(x,2,5)/dunif(x,0,1)), 0, 1)$value
D_cm <- integrate(function(y) (dbeta(y/100,2,5)/100) *
                              log((dbeta(y/100,2,5)/100)/(1/100)), 0, 100)$value
c(relative_m = D_m, relative_cm = D_cm,
  differential_m  = Hnum(function(x) dbeta(x,2,5), 0, 1),
  differential_cm = Hnum(function(y) dbeta(y/100,2,5)/100, 0, 100)) |> round(5)

##      relative_m     relative_cm  differential_m differential_cm 
##         0.48453         0.48453        -0.48453         4.12064

The relative entropy is identical both times; the differential entropy differs by \(\ln 100\). Strictly speaking, then, what one maximizes is always the relative entropy; every “maximize \(H\)” earlier in this document is shorthand for the case where the reference measure is taken to be constant (Lebesgue). That sentence is the entire content of the generalization in the next section.

7.2.3 A glossary (several names for one thing)

There is a genuine naming muddle here, worth clearing up on its own:

Term	What it refers to
reference measure / background measure / base measure	all the same thing: the \(m\) in the formula—the ruler you picked
Lebesgue measure	one specific choice of that ruler: count by ordinary length (\(m=\) const, the \(\mathrm{d}x\) of an integral)
counting measure	another specific choice: on discrete points, one notch per point (\(m(y)=1\))
density of states \(g(E)\) / degeneracy	physics’ name for the same object (the \(g\) of chapter 6)

In short: “reference measure” is the role; “Lebesgue measure” is one of the actors.

7.2.4 Three rulers worth knowing

Ruler	How it counts	What “knowing nothing” gives
equal-length (Lebesgue, \(\mathrm{d}x\))	every equal stretch counts equally—an ordinary straightedge	the uniform distribution on a bounded interval
counting (\(m(y)=1\))	one notch per integer	uniform over integers; with a fixed mean, the geometric distribution
logarithmic (\(\mathrm{d}x/x\))	every order of magnitude equally wide—like earthquake magnitude, pH, decibels	\(p\propto 1/x\), a power law

The third is worth handling directly. \(\mathrm{d}x/x\) is “equal length on a logarithmic scale” (since \(\mathrm{d}\ln x = \mathrm{d}x/x\)). Being uniform on that ruler means every decade is equally likely:

Z <- log(1000)                                # support [1, 1000], p(x) = 1/(x*Z)
sapply(1:3, function(d)
  integrate(function(x) 1/(x*Z), 10^(d-1), 10^d)$value) |>
  setNames(c("1-10", "10-100", "100-1000")) |> round(4)

##     1-10   10-100 100-1000 
##   0.3333   0.3333   0.3333

Each of the three decades holds \(1/3\)—landing between 1 and 10 is exactly as likely as landing between 100 and 1000. That is what “knowing nothing on a logarithmic scale” looks like, and on an ordinary straightedge it is the power law \(1/x\). Why does this ruler suit power laws so naturally? Because \(\mathrm{d}x/x\) is the only ruler left unchanged by \(x\to cx\) (an overall rescaling)—its unit is the factor, not the amount—and the power law is precisely the distribution with no characteristic scale.

That is all the preparation needed. The next section brings \(m\) formally into the formula and watches it light up the whole zoo at once.

7.3 The full exponential family: one template that lights up the whole zoo

Every distribution above was written as \(p\propto e^{-\lambda f(x)}\). But this is not yet the most general form—it quietly assumes one thing. Supply it and you get the object that unifies everything in statistics: the exponential family.

The hidden reference measure. Writing \(H=-\int p\ln p\,dx\) actually defaults to a reference: \(m(x)=\) constant (Lebesgue measure, “every equal-length \(dx\) counts equally”). More generally, what should be maximized is the relative entropy \(-\int p\ln\frac{p}{m}\,dx\); taking the pointwise variation with respect to \(p\) (\(\frac{\delta}{\delta p}\bigl[-p\ln\frac{p}{m}\bigr] = -\ln\frac{p}{m}-1\), setting it plus the multiplier terms to zero) solves to

\[\boxed{\,p(x)\;\propto\;m(x)\,e^{-\sum_i \lambda_i f_i(x)}\,}\]

This is the general form of the exponential family. Two parts generate the whole zoo of distributions: the reference measure \(m\) (on what “background geometry” you count) + the constraint functions \(\{f_i\}\) (which quantities you know). The preceding chapters were all special cases with \(m=\) constant and \(f\) taken as \(x\)/\(x^2\)/\(\ln x\)/\(E\).

Change the background measure and out comes a whole family of distributions.

Discrete counts \(y=0,1,2,\dots\): the most basic case—nonnegative integers + fixed mean, ordinary maximum entropy gives not Poisson but Geometric \(p_y=(1-\theta)\theta^y\) (the exponential distribution’s twin on the integer lattice). To get Poisson you must change the reference measure to \(m(y)=1/y!\); Binomial needs \(m(y)=\binom{n}{y}\). These combinatorial factors are the “background structure of counting”; the constraint cannot produce them, only \(m\) can:

Distribution	Base measure \(m(y)\)	Constraint	Result
Geometric	\(1\)	\(E[Y]=\mu\)	\((1-\theta)\theta^y\)
Poisson	\(1/y!\)	\(E[Y]=\mu\)	\(e^{-\mu}\mu^y/y!\)
Binomial	\(\binom{n}{y}\)	\(E[Y]=np\)	\(\binom{n}{y}p^y(1-p)^{n-y}\)

Scale-invariant measure \(dx/x\) (uniform in \(\ln x\), since \(d\ln x=dx/x\)): add no constraint at all and the solution is \(p\propto 1/x\)—a power law. So the power law is the distribution of “knowing nothing on a log scale,” just as the uniform is “knowing nothing on a linear scale.” The base measure is not a patch but the hidden coordinate system maximum entropy has been using all along: change it and what “most ignorant” looks like changes with it.

The big picture. \(m(x)\,e^{-\sum\lambda_i f_i}\) is not a trick of this article; it is the central object of modern statistics—the exponential family: \(\{f_i\}\) are the sufficient statistics, \(\lambda_i\) the natural parameters; generalized linear models (GLMs), conjugate priors, and the whole suite of nice maximum-likelihood properties are all built on this form. Looking back, the exponential, power-law, normal, Boltzmann, geometric, Poisson… that this article lit up one by one are just different corners of the same parent—the exponential family.

7.4 Why does nature always give these distributions?

Not because nature “chooses” these functional forms, but because:

Under constraints, the macrostate with the most microscopic realizations is the maximum-entropy state.

Nature seeks neither the lowest energy nor the highest energy, but the most easily realized state—the distribution with the most microscopic arrangements. The exponential decay \(e^{-\lambda f(x)}\) recurs because it is the general solution of the constrained optimization problem.

8 Overall summary

One main thread runs through the whole article: optimize a quantity + impose some constraints → a uniquely determined structure.

1. What is optimized: entropy \(H = -\int p\ln p\,dx\). It is not arbitrarily chosen—the counting argument proved \(\ln W \approx N\cdot H\), and the maximum-entropy distribution is the one with overwhelmingly the most microscopic arrangements and therefore all but certain to appear (\(W \approx e^{NH}\)).
2. How to optimize: Lagrange multipliers / calculus of variations. Absorb the constraints into the Lagrangian, differentiate pointwise with respect to the distribution, and the general solution drops out in one step: \(p(x) \propto e^{-\lambda f(x)}\).
3. What the constraints are: determined by your known information (measurement or conservation law). Change the constraint function \(f(x)\) and you change the distribution:

Constraint \(f(x)\)	Resulting distribution	Physical identity of the multiplier	Example
none (bounded only)	uniform distribution	—	—
\(x\) (mean, nonnegative)	exponential distribution	—	isothermal atmosphere, radioactive decay
\(\ln x\) (geometric mean, \(x\ge x_{\min}\))	power-law / Pareto	—	city sizes, wealth tails, network degrees
\(x,\ x^2\) (mean + variance)	normal distribution	—	chemical-bond vibration, Brownian diffusion
\(E(x)\) (energy)	Boltzmann distribution	\(\lambda = 1/k_BT\) (temperature)	Maxwell speeds, Arrhenius

Exponential, power-law, normal, Boltzmann—these four seemingly unrelated distributions are therefore four faces of the same equation under different constraints (the uniform being the “no constraint” baseline). And this “optimization + constraints” template goes far beyond this: the Lagrange multiplier is a constraint force in mechanics, temperature/chemical potential/pressure in thermodynamics, and the dual variable in optimization—it is a common axis connecting mathematics, physics, and chemistry.

One line to take away: the exponential family is not nature’s aesthetic preference but the mathematical consequence of one inferential situation—a fixed reference measure, and knowledge of only a handful of average quantities.

The weight of that sentence has to be placed carefully: its necessity belongs to inference, not to nature. Loosen any of the three premises (fixed reference measure, finitely many known constraints, maximize entropy) and the conclusion fails—and nature does contain distributions outside the exponential family: Student’s \(t\), Cauchy, Lévy stable laws, mixtures, and families whose support moves with the parameters (Pareto with unknown \(x_{\min}\), Weibull with unknown shape). The converse matters just as much: distributions that are exponential families frequently arise for reasons having nothing to do with honest inference—the normal grows out of the central limit theorem (additive aggregation), power laws out of preferential growth, Kesten processes, or critical branching.

So one and the same distribution often has two independent provenances: maximum entropy answers “what should I assume given this much information,” while mechanism answers “what will it actually grow into.” The two routes coincide remarkably often (the normal, the exponential and the power law each have both kinds of explanation), yet they are answers to different questions—and that coincidence is the part most worth turning over.

---

By Dr. Yang Hu. English translation drafted with Claude (Opus 4.8), reviewed by the author.

First published: 2026-07-12. Last updated: 2026-08-26 15:37 EDT.

References

Beggs, John M., and Dietmar Plenz. 2003. “Neuronal Avalanches in Neocortical Circuits.” The Journal of Neuroscience 23 (35): 11167–77. https://doi.org/10.1523/JNEUROSCI.23-35-11167.2003.

Boltzmann, Ludwig. 1877. “Ueber Die Beziehung Zwischen Dem Zweiten Hauptsatze Der Mechanischen w"armetheorie Und Der Wahrscheinlichkeitsrechnung Respektive Den s"atzen "uber Das w"armegleichgewicht.” Sitzungsberichte Der Kaiserlichen Akademie Der Wissenschaften. Mathematisch-Naturwissenschaftliche Classe 76: 373–435. https://doi.org/10.3390/e17041971.

Broido, Anna D., and Aaron Clauset. 2019. “Scale-Free Networks Are Rare.” Nature Communications 10: 1017. https://doi.org/10.1038/s41467-019-08746-5.

Cercignani, Carlo. 2006. Ludwig Boltzmann: The Man Who Trusted Atoms. Oxford University Press.

Clauset, Aaron, Cosma Rohilla Shalizi, and Mark E. J. Newman. 2009. “Power-Law Distributions in Empirical Data.” SIAM Review 51 (4): 661–703. https://doi.org/10.1137/070710111.

Gabaix, Xavier. 2009. “Power Laws in Economics and Finance.” Annual Review of Economics 1: 255–94. https://doi.org/10.1146/annurev.economics.050708.142940.

Harris, Theodore E. 1963. The Theory of Branching Processes. Vol. 119. Grundlehren Der Mathematischen Wissenschaften. Springer. https://doi.org/10.1007/978-3-642-51850-3.

Jaynes, Edwin T. 1957. “Information Theory and Statistical Mechanics.” Physical Review 106 (4): 620–30. https://doi.org/10.1103/PhysRev.106.620.

Kesten, Harry. 1973. “Random Difference Equations and Renewal Theory for Products of Random Matrices.” Acta Mathematica 131: 207–48. https://doi.org/10.1007/BF02392040.

Moivre, Abraham de. 1733. Approximatio Ad Summam Terminorum Binomii (a+b)^n in Seriem Expansi. https://doi.org/10.5281/zenodo.18913559.

Planck, Max. 1901. “Ueber Irreversible Strahlungsvorg"ange.” Annalen Der Physik 311 (12): 818–31. https://doi.org/10.1002/andp.19013111210.

Shannon, Claude E. 1948. “A Mathematical Theory of Communication.” Bell System Technical Journal 27 (3–4): 379–423, 623–56. https://doi.org/10.1002/j.1538-7305.1948.tb01338.x.

Stirling, James. 1730. Methodus Differentialis: Sive Tractatus de Summatione Et Interpolatione Serierum Infinitarum. G. Strahan. https://books.google.com/books?id=71ZHAAAAYAAJ.

Tribus, Myron, and Edward C. McIrvine. 1971. “Energy and Information.” Scientific American 225 (3): 179–88. https://doi.org/10.1038/scientificamerican0971-179.

West, Geoffrey B., James H. Brown, and Brian J. Enquist. 1997. “A General Model for the Origin of Allometric Scaling Laws in Biology.” Science 276 (5309): 122–26. https://doi.org/10.1126/science.276.5309.122.

---

Comments & discussion

Anything unclear, anything wrong, or an example worth adding — please leave a note below (a GitHub account is all it takes). Comments support Markdown, code blocks and $LaTeX$.