从熵到分布：指数、幂律、正态与玻尔兹曼的统一逻辑
胡杨博士
首次发布：2026-07-12 | 最后更新：2026-08-26 16:24 EDT
本文的核心论点：指数分布、幂律分布、正态分布、玻尔兹曼分布不是四个独立的分布，而是同一个原理（最大熵）在不同约束下的四种表现。
\[\text{最大化 } H = -\!\int p(x)\ln p(x)\,dx \quad \text{subject to constraints} \quad \Longrightarrow \quad p(x) \propto e^{-\lambda\, f(x)}\]
把四个 \(f(x)\) 代进这一个式子，四个分布就依次掉出来：
| 分布 | 约束的量 \(f(x)\) | 解析形式（\(e^{-\lambda f(x)}\) \(\xrightarrow{\text{归一化}}\) \(p(x)\)） |
|---|---|---|
| 指数分布 | \(x\) —— 均值（\(x \ge 0\)） | \(e^{-\lambda x} \;\to\; \lambda e^{-\lambda x}\) |
| 幂律 / Pareto | \(\ln x\) —— 几何均值（\(x \ge x_{\min}\)） | \(e^{-\lambda \ln x} = x^{-\lambda} \;\to\; (\lambda-1)\,x_{\min}^{\lambda-1}\,x^{-\lambda}\) |
| 正态分布 | \(x,\ x^2\) —— 均值与方差 | \(e^{-\lambda_1 x-\lambda_2 x^2} \;\to\; \dfrac{1}{\sqrt{2\pi}\,\sigma}\,e^{-(x-\mu)^2/2\sigma^2}\) |
| 玻尔兹曼分布 | \(E(x)\) —— 能量 | \(e^{-\lambda E(x)} \;\to\; \dfrac{1}{Z}\,e^{-E(x)/k_BT}\) |
第二行值得盯一眼：\(e^{-\lambda\ln x}\) 就是 \(x^{-\lambda}\)——把约束从”\(x\) 的均值”换成”\(\ln x\) 的均值”，指数分布当场变成幂律。同理第三行多加一个 \(x^2\) 就配出高斯的平方项，第四行把 \(\lambda\) 认成 \(1/k_BT\) 就是物理学家的玻尔兹曼因子。四个分布之间没有跨越，只有 \(f(x)\) 换了一次。
1 为什么是最大熵？
在最大化熵之前，先回答一个更根本的问题：凭什么要最大化熵？ 有两条独立的路都通向它——物理的计数论证，和统计的无偏推断论证。
1.1 计数论证：熵 = 排列数的对数
做 \(N\) 次观测，每次结果落入 \(k\) 个类别之一。记 \(n_i\) 为类别 \(i\) 出现的次数（\(i = 1, \dots, k\)，满足 \(n_1 + n_2 + \cdots + n_k = N\)）。宏观状态是计数向量 \((n_1, \dots, n_k)\)——只记”每类出现多少次”；微观状态是具体的观测序列——还记”哪一次是哪一类”。一个宏观状态对应的微观排列数是多项式系数：
\[W = \frac{N!}{n_1!\,n_2!\cdots n_k!}\]
抛硬币是 \(k=2\) 的特例：类别只有正/反，\(n_1 = k_{\text{正}}\)，\(n_2 = N - k_{\text{正}}\)，于是
\[W = \frac{N!}{k_{\text{正}}!\,(N-k_{\text{正}})!} = \binom{N}{k_{\text{正}}}\]
证明 \(\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{，两个线性项相消} \\[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{，常数 } \ln N \text{ 乘进求和} \\[2pt] &= \sum_i n_i \bigl(\ln N - \ln n_i\bigr) && \text{同一求和指标，合并并提出 } n_i \\[2pt] &= -\sum_i n_i \ln\frac{n_i}{N} && \textstyle\ln N - \ln n_i = -\ln\frac{n_i}{N} \text{（对数商法则），负号提出} \\[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} \]
即 \(W \approx e^{N H}\)：熵是每次观测平摊到的”排列数指数”。两个分布的熵哪怕只差一点点，排列数就差 \(e^{N\cdot\Delta H}\) 倍——\(N\) 大时是天文数字。所以”熵最大的分布”不是抽象偏好，而是压倒性地更可能出现的宏观状态。
## ── 计数直觉：熵大 = 微观排列数多 = 压倒性地更可能 ──────────
par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))
# 左图：抛 30 次硬币，正面次数 k 对应的排列数 W = C(30, k)
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)
# 右图：(1/N) ln C(N, pN) 随 N 增大收敛到二元熵 H(p)
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)
左图：每根柱子是一个宏观状态（正面占比 \(k/N\)），柱高是它对应的微观排列数 \(W = \binom{30}{k}\)。注意 \(k/N = 0\) 和 \(1\) 的柱子画了但看不见——高度是 \(W = 1\)（全反/全正只有一种排法），在 1.55 亿的纵轴下肉眼就是零。这正是要传达的信息：\(k/N = 1/2\) 的排列数是两端的 1.5 亿倍，公平硬币下你”看到”一半正面，不是有什么力量在推它，纯粹是排列数碾压。右图：\(\frac{1}{N}\ln W\) 随 \(N\) 增大收敛到 \(H(p)\)，正是上面 Stirling 证明的数值验证。
1.2 峰有多宽？二阶展开与 \(1/\sqrt{N}\) 窗口
上图有个容易起疑的地方：说好的”碾压”，中间几根柱子怎么差得不大？比如 \(\binom{30}{14}/\binom{30}{15} = 0.94\)，只差 6%。这不是矛盾，是二阶展开的必然。
\(W_1/W_2 = e^{N\Delta H}\) 说的确实是 \(W\) 本身的比值，但关键在 \(\Delta H\) 有多大。\(H(p)\) 在 \(p = 1/2\) 处是光滑极大值，一阶导数为零，Taylor 展开从二阶开始（记 \(\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\]
于是峰附近：
\[\frac{W(p)}{W_{\max}} \approx e^{-2N\delta^2}\]
两个结论直接读出来：
- 相邻柱子差不大：\(\delta = 1/30\) 时 \(e^{-2N\delta^2} \approx e^{-1/15} \approx 0.94\)——你看到的”中间平”就是一阶导为零的表现；
- 窗口宽度 \(\sim 1/\sqrt{N}\)：\(e^{-2N\delta^2}\) 掉到可忽略需要 \(\delta \gtrsim 1/\sqrt{N}\)。这正是 CLT 的 \(\sqrt{N}\) 涨落——排列数计数和中心极限定理在这里是同一件事。
所以”碾压”是对固定的宏观差异（固定 \(\Delta H\)）、随 \(N\) 增大而言的；\(N\) 固定时峰附近永远有一个 \(O(1/\sqrt{N})\) 宽的”没被碾压”窗口，\(N\to\infty\) 时窗口本身收缩成零。
## ── 碾压是 N 的函数：窗口收缩 + 固定 ΔH 的指数塌缩 ──────────
par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))
# 左图：W/W_max 随 N 增大收缩成尖峰（宽度 ~ 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")
# 右图：固定宏观差异 p=0.6 vs p=0.5，W 比值随 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)
左图：把每个 \(N\) 的排列数都除以自己的峰值，\(N = 30 \to 3000\)，窗口按 \(1/\sqrt{N}\) 收缩。右图：固定比较 \(p = 0.6\) 对 \(p = 0.5\)（\(\Delta H \approx 0.0201\)），\(W\) 比值在对数轴上是一条直线 \(e^{-N\Delta H}\)——\(N=30\) 时只差一半（0.56，毫不起眼），\(N=1000\) 时已是 \(10^{-9}\)。精确二项计算（蓝）和 \(e^{-N\Delta H}\)（红虚线）几乎重合。
1.3 历史注记：同一个 \(H\)，两次被发现
上面的推导看着像信息论，历史上却是物理先走通的——\(-\sum p_i \ln p_i\) 这个量被独立发现了两次，路径完全不同：
1730–1733 年，de Moivre 与 Stirling：de Moivre 在二项分布近似中得到 \(n! \sim C\,n^{n+1/2}e^{-n}\) 的结构；Stirling 在 Methodus Differentialis（1730）确定 \(C=\sqrt{2\pi}\)。de Moivre 随后在 1733 年的小册子里把它用于二项分布的正态近似。工具就位。(Moivre 1733; Stirling 1730)
1877 年，Boltzmann：本节的计数思路直接承接他这年论文——把气体分子的能量切成离散小份，数每种分配的排列数（他称 Komplexionen，即 \(W\)），再用 Stirling 分析最可能的宏观状态。\(S = k_B \ln W\) 后来刻在他的墓碑上；除以 \(N\)，就出现 \(-\sum p_i \ln p_i\) 的形式。(Boltzmann 1877)
数学在今天看是本科习题，但公式的难度 ≠ 洞见的难度——每一步在当时都带有争议：统计解释的第二定律引出 Loschmidt 可逆性与 Zermelo 回归的反驳；原子实在性又遭 Mach、Ostwald 等人质疑。Boltzmann 1906 年自杀，未能亲见 Perrin 随后以布朗运动实验显著推动原子论的接受。Planck 在黑体辐射工作中把离散能量元保留下来，给出了后来标准写法中的常数 \(k\) 并以 Boltzmann 命名。(Cercignani 2006; Planck 1901)
1948 年，Shannon（《A Mathematical Theory of Communication》）：走公理化路线——要求不确定性度量对 \(p\) 连续、等概率时随类别数单调增、分步选择时可加权分解，证明唯一解是 \(-K\sum p_i \log p_i\)。推导里没有排列数、没有 Stirling；但典型序列又绕回计数：长度 \(N\) 的序列里”典型”的约有 \(e^{NH}\) 条。(Shannon 1948)
命名轶事：后来广为流传的版本称，von Neumann 建议 Shannon 把这个量叫”熵”，因为”没人真正懂熵是什么”。这段话可追溯到 1971 年的转述，而非可核验的谈话记录；把它当作好玩的二手轶事即可。(Tribus and McIrvine 1971)
1957 年，Jaynes（《Information Theory and Statistical Mechanics》）：把两条路正式焊在一起——统计力学的最大熵不是物理定律，而是推断原理的特例。这就是下一小节。(Jaynes 1957)
1.4 Jaynes：最大熵 = 最诚实的推断
第二条路不谈物理，谈推断。(Jaynes 1957) 当你只知道一个分布的部分信息（比如只知道均值），却要选出一个完整的分布时：
- 选熵最大的那个 = 只用上已知的约束，不夹带任何额外假设；
- 选其他任何分布 = 偷偷宣称你知道一些其实不知道的信息（比如”分布是双峰的”）。
所以最大熵不是自然界的神秘偏好，而是约束之外保持最大无知的唯一一致做法。物理系统”选”它是因为排列数碾压，统计学家”选”它是因为诚实——两条路殊途同归，都指向同一个变分问题：
\[\max_p \; H[p] \quad \text{s.t. 已知约束} \;\Longrightarrow\; p(x) \propto e^{-\lambda f(x)}\]
后面的第3章到第6章就是把不同的 \(f(x)\) 代进去。
自行观看（拓展）。 Reinventing Entropy | Compression is Intelligence Part 1（3Blue1Brown）从压缩的角度重新搭建熵的直觉；它比本文的最大熵主线更宽，可在读完本节后观看。
2 熵：为什么是这些分布，而不是别的？
2.1 信息熵的定义
对连续分布 \(p(x)\)，微分熵定义为：
\[H[p] = -\int p(x) \ln p(x)\, dx\]
\(H\) 越大，分布越”分散”、越”不确定”、越”没有偏见”。
（这个定义里其实藏了一把”尺子”：所谓”分散”是相对什么在量?连续情形下换个单位，\(H\) 的数值就会变。这件事到末章《先说清那把尺子》才需要正面处理，在那之前按常规读法即可。）
记号注记：为什么信息熵写 \(H\) 而不是 \(S\)。 这个字母是 Shannon 从 Boltzmann 的 H 定理借来的——他在 1948 年的论文里明说了，\(H\) 就是 H 定理里那个 \(H = \int f\ln f\,dv\)（Boltzmann 证明它单调减，所以 \(S = -k_B H\)；Shannon 借了字母、把符号翻了过来）。至于 Boltzmann 为何用 H，一说它本是希腊大写 Eta（Η），因为熵在部分传统里记作小写 \(\eta\)，而大写 Eta 在印刷上与拉丁 H 无从分辨。
要紧的是两者不是同一个量，差着一个 \(k_B\) 和量纲：
| 记号 | 定义 | 量纲 | |
|---|---|---|---|
| 信息熵（本文主用） | \(H\) | \(-\sum p\ln p\) | 无量纲（nat；换 \(\log_2\) 则为 bit） |
| 热力学熵 | \(S\) | \(-k_B\sum p\ln p\) | J/K |
\[S = k_B H\]
本文按领域惯例分用两个字母：讲推断时写 \(H\)，讲热力学（温度、玻尔兹曼分布）时写 \(S\)，两处用上式换算。物理学界内部也常直接用 \(S\) 写这同一个泛函（Gibbs 熵、Tsallis 熵 \(S_q\)），所以 \(H\) 与 \(S\) 是惯例之差，不是对错之分。另外提醒一句：热化学里的 \(H\) 是焓，与熵无关，本文完全不涉及。
2.2 读懂被积函数：\(p(x) \times (-\ln p(x))\)
定义里的乘积 \(-p\ln p\) 值得拆开看。改写成 \(H = \int p(x)\cdot\bigl(-\ln p(x)\bigr)\,dx = E\bigl[-\ln p(X)\bigr]\)：
- \(-\ln p(x)\) 是在 \(x\) 处的信息量（惊讶度）：\(p\) 越小越惊讶，尾部趋于 \(+\infty\)；
- \(p(x)\) 是概率权重：这个惊讶度实际发生的频率；
- 熵 = 平均惊讶度。
## ── 拆开被积函数：p(x) × (-ln p(x)) ──────────────────────────
x_e <- seq(-4, 4, length.out = 2000)
p_e <- dnorm(x_e) # 概率权重
s_e <- -log(p_e) # 信息量（惊讶度）
g_e <- p_e * s_e # 熵密度 = 两者乘积
pts <- c(0, 1.5, 3) # 三个跟踪点：中心 / 中段 / 尾部
pt_cols <- c("tomato", "purple", "forestgreen")
par(mfrow = c(1, 3), mar = c(4.6, 4.8, 3.5, 1))
# 面板 1：概率权重 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)
# 面板 2：信息量 -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)
# 面板 3：乘积 = 熵密度，阴影面积 = 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)
跟踪三个点（标准正态）：
- 红点 \(x=0\)：权重最大（0.399），但惊讶度最小（0.92）——乘积中等偏大；
- 紫点 \(x=1.5\)：权重和惊讶度都中等——乘积仍然可观；
- 绿点 \(x=3\)：惊讶度很大（5.4），但权重几乎为零（0.004）——乘积 \(\approx 0\)。
尾部不失控是因为 \(\lim_{p\to 0} p\ln p = 0\)：\(p\) 趋于 0 的速度（指数级）永远快过 \(\ln p\) 发散的速度（对数级）。所以熵的贡献主要来自”中等概率”区域，第三幅图阴影面积就是 \(H = \tfrac{1}{2}\ln(2\pi e) \approx 1.419\)。
2.3 熵随分布展宽而增大
## ── 直觉：熵衡量分布的"展开程度" ──────────────────────────
# 辅助函数：数值计算微分熵
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]
# 不同 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))
# 左图：分布越宽 → 熵越大
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)
# 右图：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)
直觉总结：
- 把概率集中在一个点 → \(H \to -\infty\)（完全确定，熵最低）
- 均匀展开 → \(H\) 最大（完全不确定，熵最高）
- 在有约束的情况下，分布不能完全均匀，但可以在约束允许的范围内尽量”展开”
2.4 最大熵原理
当你只知道关于数据的某些统计量（均值、方差……），最诚实的做法是选择满足这些约束、但其他方面尽可能”不做假设”的分布——即熵最大的分布。
（英文 Maximum Entropy，缩写 MaxEnt，本文图内英文标注用它。）
这不是一个任意的审美偏好，而是一个逻辑必然：
- 如果你选了一个熵更低的分布，你就在暗中假设你实际上不知道的信息
- 最大熵分布 = 在约束之外不注入任何额外信息的分布
2.5 拉格朗日乘子：零基础版
本节和下一节（拉格朗日乘子的几何与推导）是纯方法铺垫。如果你已经学过拉格朗日乘子法，或暂时不想看推导，可以放心跳过这两节，直接到《约束决定分布》看结论，甚至直接进入后面的具体分布（从《指数分布》那章开始）——不影响后面的阅读。
下一小节要用拉格朗日乘子法（Lagrange multipliers），先把这个工具本身讲清楚。核心就一张图，数学只有一行。
问题形态。普通极值问题：\(f(x,y)\) 哪里最大？——导数为零，\(\nabla f = 0\)。带约束的极值问题：在满足 \(g(x,y) = c\) 的前提下 \(f\) 哪里最大？你只能在约束曲线上走，最优点处 \(f\) 的导数一般不是零——只是你被约束拦住了。
几何直觉（整个方法的灵魂）。沿约束曲线走，盯着 \(f\)：
- 若 \(\nabla f\)（\(f\) 增长最快的方向）在曲线切向上有分量，沿曲线挪一步就能改进——还没到最优；
- 到了最优点，\(\nabla f\) 必须垂直于约束曲线（切向分量为零，无路可走）。
而 \(\nabla g\) 永远垂直于曲线 \(g = c\)（沿线 \(g\) 不变）。两个向量垂直于同一条曲线 → 平行：
\[\boxed{\nabla f = \lambda\, \nabla g}\]
\(\lambda\) 就是拉格朗日乘子——它的全部含义就是”最优点处两个梯度平行”的比例系数。
操作配方。构造拉格朗日函数把约束吃进目标函数：
\[L(x, y, \lambda) = f(x,y) - \lambda\,\bigl(g(x,y) - c\bigr)\]
对所有变量（含 \(\lambda\)）求偏导设零：前两个方程给出 \(\nabla f = \lambda\nabla g\)，\(\partial L/\partial\lambda = 0\) 自动还原约束。带约束问题变成了无约束问题——这就是它好用的原因。
例：周长 20 的矩形，面积最大？最大化 \(f = xy\)，约束 \(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\]
正方形，面积 25。
下图把这件事画成一座山：曲面是完整的 \(f = xy\)，整座山一直朝 \((10,10)\) 方向涨上去。蓝色环线是等高线（地形图画法，每条线上 \(f\) 相同）；\(\nabla f\) 的方向由此一眼可读——垂直穿过等高线、指向更高一环，这就是最陡上坡方向。约束线抬升到曲面上是横穿山坡的一条路（红色），约束优化 = 只许沿这条路走，找路上的最高点。关键看紫色点 \((5,5)\)：\(\nabla f\) 在这里并不为零——绿色箭头垂直于等高线指向上坡，\(f\) 还能继续变大，但那个方向离开了约束线，不许走；为零的只是 \(\nabla f\) 沿路的切向分量（右图：路的高度剖面在此处坡度为零）。而在紫色点处约束线恰好与等高线 \(f = 25\) 相切——路贴着等高线走的瞬间，\(f\) 不增不减，正是极值。
图中还画了两样东西。半透明的竖墙是约束 \(g = 10\) 在三维里的形象：\(g(x,y) = x+y\) 不含 \(z\)，所以 \(\{x + y = 10\}\) 在三维里是一面 \(z\) 自由伸展的墙（这里只画到山体表面为止），红路就是墙切进山体的截口；地面虚线是墙脚，即约束集本身。橙色箭头是 \(\nabla g\) 升维成 \((1,1,0)\) 后的样子：它是这面墙的唯一法向——垂直于三维里的一条曲线有无穷多个方向（一整个法平面），但垂直于一张曲面只有一个方向，所以”\(\nabla g\) 垂直于约束”用墙来想最不容易误会。地面二维版是同一句话：\(\nabla g = (1,1)\) 只有两个分量、出不了地面，在平面内垂直于虚线，方向同样唯一。
梯度是几维的？ 分量数 = 输入个数。\(f(x,y) = xy\) 两个输入，所以 \(\nabla f = (y,\ x)\) 是二维向量，\(\nabla f(5,5) = (5,\ 5)\)；\(z = 25\) 是输出，不占分量。图像点 \((5,5,25)\) 在三维，但 \(f\) 和它的梯度都住在二维定义域（“地图”）里。方向 \((1,1)\) 给出最陡上坡，模长给出坡度：\(|\nabla f| = \sqrt{50} \approx 7.07\)，即沿最陡方向水平走 1，\(f\) 涨 7.07。
亲手验证一遍。沿 \((1,1)\) 方向走，先归一化成单位向量 \(\bigl(\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2}\bigr)\)（向量 \((1,1)\) 长 \(\sqrt2\)，不归一化 \(t\) 就不是走过的距离），走过距离 \(t\) 后位置是 \(x(t) = y(t) = 5 + \tfrac{t}{\sqrt2}\)，于是：
\[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\]
同一结果也可由方向导数公式一步得到：\(\nabla f \cdot \hat u = (5,5)\cdot\bigl(\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2}\bigr) = \tfrac{10}{\sqrt2}\)——\(\hat u\) 与 \(\nabla f\) 同向时点积最大，恰为 \(|\nabla f|\)。“梯度模长 = 最陡坡度”不是定义，是算出来的。
据此校准图例：绿箭头是贴着山坡的可视化向量，其水平投影才是 \(\nabla f\)；\(\nabla g\) 画成 \((1,1,0)\) 同理（补 0 入图）。真正三维的梯度属于另一个函数——\(F(x,y,z) = z - xy\) 的 \(\nabla F = (-y, -x, 1)\)，垂直于山坡曲面本身。一般规律：\(n\) 元函数的梯度是 \(n\) 维向量，垂直于 \((n-1)\) 维等值集。拉格朗日条件用的是 \(\nabla f\)。
## ── Lagrange 乘子的几何：约束线是曲面上的一条"山路" ──────────
par(mfrow = c(1, 2), mar = c(1, 1, 3.5, 0.5))
# 左图：曲面 z = xy，约束 x + y = 10 抬到曲面上是一条山路
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)
# 补全边框立方体：persp 只画了靠近原点的隐藏边，把 x=10 / y=10 侧补齐
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) # (10,10) 角的竖边
bx(10, 0, 0, 10, 10, 0) # 底面后边 x = 10
bx( 0, 10, 0, 10, 10, 0) # 底面后边 y = 10
bx(10, 0, 100, 10, 10, 100) # 顶面后边 x = 10
bx( 0, 10, 100, 10, 10, 100) # 顶面后边 y = 10
# 竖墙：g = 10 沿 z 方向延伸（x+y=10, z 自由），墙切山的截口 = 红路
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)：墙的唯一法向（水平，垂直于墙面）
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)
# 山体表面的等高线（地形图环线）：grad f 处处垂直于它们
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)
# 约束线（地面上）
lines(trans3d(x_p, 10 - x_p, 0, pm), col = "gray30", lwd = 2, lty = 2)
# 山路：约束线抬升到曲面上，高度 f = x(10-x)
lines(trans3d(x_p, 10 - x_p, x_p * (10 - x_p), pm), col = "tomato", lwd = 3.5)
# 非最优点 (2,8) 与最优点 (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)
# 最优点垂线
lines(trans3d(c(5, 5), c(5, 5), c(0, 25), pm), col = "purple", lty = 3, lwd = 1.5)
# grad f 在 (5,5) 不为零：指向上坡（离开约束线）
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)
# 右图：沿约束线走，f 的值 = x(10-x)（山路的侧面展开）
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)
\(\lambda\) 的含义：约束的”价格”。一个漂亮的定理：最优值对约束的敏感度就是乘子本身，
\[\frac{df^*}{dc} = \lambda\]
上例 \(\lambda = 5\)：周长的一半从 10 放松到 11，最大面积涨约 5（精确值 \(5.5^2 - 5^2 = 5.25\)）。（\(\lambda\) 的数值绑定在约束的写法上：把 \(g\) 乘 2，\(\lambda\) 减半，\(df^*/dc\) 同步减半，定理在任何写法下自洽；量纲是 \([f]/[g]\)。）
这个”约束的价格”在各领域都有自己的名字：
| 领域 | \(\lambda\) 的名字 | 例子 |
|---|---|---|
| 优化/数学 | 对偶变量、KKT 乘子 | 线性规划对偶问题的解 |
| 分析力学 | 约束力 | 钢丝给珠子的法向支持力、绳的张力 |
| 统计物理 | 共轭强度量 | \(\beta = 1/k_BT\)（能量约束） |
| 化学 | 化学势 \(\mu\) | 粒子数约束 |
| 经济学 | 影子价格 | 预算约束 = 边际效用 |
力学那行最直观：小球被约束在斜面上，拉格朗日方程解出的 \(\lambda\) 就是斜面顶住小球的法向力——“支持力垂直于接触面”和”\(\nabla g\) 垂直于约束集”是同一句话。热力学的强度量–广延量配对表（\(T\)–\(S\)、\(p\)–\(V\)、\(\mu\)–\(N\)）本质上是一张拉格朗日乘子登记表：每个强度量都是某个守恒量约束的乘子，这也是它们天生不随系统变大的原因——乘子是比率 \(\partial f^*/\partial c\)。《玻尔兹曼分布》那章会看到其中第一行：温度就是能量约束的影子价格，\(\partial S/\partial E = 1/T\)。
推广到分布：最大熵问题里”变量”不是两个数，而是每一点的密度值 \(p(x)\)——无穷多个变量，每个 \(x\) 一个。配方不变（这一步叫变分法），下一小节就做这件事。物理里这套方法无处不在，因为物理问题几乎全是”某量取极值 + 守恒律当约束”，而乘子往往就是有名字的物理量：温度（能量约束）、化学势（粒子数约束）、压强（体积约束）。
2.6 拉格朗日乘子推导：分布的一般形式
最大化 \(H = -\int p \ln p\, dx\)，约束条件：
- 归一化：\(\int p(x) dx = 1\)
- 某个函数的期望值固定：\(\int p(x) f(x) dx = \langle f \rangle\)
（约束 2 不是推导出来的，是问题的输入：你测到或守恒律给定了某个平均量，“知道一个平均量”写成数学就是期望值形式——\(f(x)\) 指明测的是什么量，\(\langle f\rangle\) 是测出来的那个数。对照：
- 测到样本均值 \(\bar x\)：取 \(f(x) = x\)，约束为 \(\int p(x)\,x\,dx = \bar x\)，即 \(\langle f\rangle = \bar x\)；
- 测到方差 \(s^2\)（连同均值 \(\bar x\)）：再加一条 \(f(x) = x^2\)，约束为 \(\int p(x)\,x^2\,dx = s^2 + \bar x^2\)（二阶矩），即 \(\langle f\rangle = s^2 + \bar x^2\)；
- 能量守恒、平均每个粒子分到 \(\bar E\)：取 \(f(x) = E(x)\)，约束为 \(\int p(x)\,E(x)\,dx = \bar E\)，即 \(\langle f\rangle = \bar E\)。
\(f\) 和 \(\langle f\rangle\) 都由已知信息决定，最大熵原理只负责已知之外不多假设——所以换一个 \(f\) 就换一个分布。）
构造拉格朗日函数：
\[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)\]
对 \(p(x)\) 求变分并令其为零：
“求变分”是什么意思？ 普通微积分优化一个数，这里优化一整条函数（找哪条 \(p\) 让 \(H[p]\) 最大）——这就是变分法。实操上只需一条规则：被积函数 \(F(p)\) 对 \(p\) 当普通变量求导即可，
\[\frac{\delta}{\delta p(x)} \int F(p)\, dx = \frac{\partial F}{\partial p}\]
用到熵项：\(F = p\ln p\)，\(\frac{\partial}{\partial p}(p\ln p) = \ln p + 1\)（乘积法则，\(+1\) 来自 \(p\cdot\frac{1}{p}\)）。所以 \(-\int p\ln p\) 贡献 \(-\ln p - 1\)，归一化项贡献 \(-\lambda_0\)，期望项贡献 \(-\lambda_1 f(x)\)，相加即下式。（\(+1\) 最后并进 \(Z = e^{1+\lambda_0}\)，不影响解。）
\[\frac{\delta L}{\delta p} = -\ln p(x) - 1 - \lambda_0 - \lambda_1 f(x) = 0\]
解出：
\[\boxed{p(x) = \frac{1}{Z}\, e^{-\lambda_1 f(x)}}\]
其中 \(Z = e^{1+\lambda_0}\) 是归一化常数（配分函数），\(\lambda_1\) 由约束 \(\langle f\rangle\) 决定。
这就是全部的推导。 所有分布都从这一个公式出来，区别只在于 \(f(x)\) 是什么。
2.7 约束决定分布
par(mfrow = c(2, 2), mar = c(4, 4, 3.5, 1))
# ── 1. 无约束（有界区间） → 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. 固定均值（非负） → 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. 固定均值 + 方差 → 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. 固定对数均值 → 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)
| 约束 | \(f(x)\) | \(\lambda\) | 最大熵分布 \(p(x)\) |
|---|---|---|---|
| 无约束（有界） | — | — | Uniform：\(\dfrac{1}{b-a}\) |
| \(E[X] = \mu\)，\(X \geq 0\) | \(x\) | \(1/\mu\) | Exponential：\(\lambda e^{-\lambda x}\) |
| \(E[\ln X]\) 固定，\(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}\) |
所有这些分布来自同一个方程的不同解。 区别只在于你”知道什么”（约束）。这里的”约束”不只包括均值、方差、能量这类期望条件，也包括支撑集本身：均匀分布要先知道变量被限制在 \([a,b]\)，指数分布要先知道 \(X\ge0\)，幂律分布要先知道 \(X\ge x_{\min}>0\)。最大熵问题从来不是”凭空找一个分布”，而是：
\[\text{给定支撑集 + 给定若干统计量} \quad\Longrightarrow\quad \text{在这些信息内熵最大}\]
3 指数分布 = 最大熵 + 非负 + 均值约束
3.1 推导
约束：\(X \geq 0\)，\(E[X] = \mu\)。这里 \(X\ge0\) 本身就是一条支撑约束：等待时间、寿命、高度、能量这类量不能为负；在这个半轴上再固定均值，才得到指数分布。代入通解：
\[p(x) = \frac{1}{Z} e^{-\lambda x} = \lambda e^{-\lambda x}, \quad \lambda = 1/\mu\]
归一化常数 \(Z = 1/\lambda\)。微分熵 \(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{拆开，}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} \]
单位是 nat（用 \(\ln\)）。均值越大分布越展开，\(H\) 随之增大；\(\mu < 1/e\) 时 \(H < 0\)——微分熵可负，与离散熵不同。
3.2 数值验证：候选分布中，指数的熵最高
在均值 \(= 2\) 的所有非负分布中比较熵：
# 在 [0, max] 网格上数值计算微分熵
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 # 固定均值
# 均值 = 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) # 均值 2 的均匀分布：熵反而更小
)
entropies <- sapply(dists, h_pos)
h_theory <- 1 + log(mu) # Exp(1/mu) 的解析熵
par(mfrow = c(1, 2), mar = c(4, 7, 3.5, 1))
# 左图：密度对比
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)
# 右图：熵的比较（横向条形图）
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)
Exp(1/2) 的红色条最长——在均值 \(= 2\) 的所有非负分布中，指数分布的熵确实最大。
Gamma(k) 随着 \(k\) 增大（形状越像正态），熵越来越低，因为分布越来越”确定”。
注意最下面那条 Uniform[0,4]（均值也是 2）：它的熵 \(\ln 4 \approx 1.386\)，比指数的 \(1.693\) 还小。这看着反直觉——“均匀 = 最大熵”不是常说的吗？那条结论只在支撑集固定有界、且无其他约束时成立（“\(X\) 一定在 \([0,4]\) 内，别的不知道”→ 均匀）。这里支撑是整个 \([0,\infty)\)、约束是均值：\(U[0,4]\) 在 \(x>4\) 处概率硬性为零，等于偷偷断言”\(X\) 绝不超过 4”——一条你并不知道的额外信息，反而扣熵。指数不设上限、铺满半轴，除均值外不假设任何东西，所以熵更高。
3.3 严格证明：固定均值下，指数熵最大（不只是”图上最高”）
上图只比了 7 条。要证的命题是：在均值固定为 \(\mu\) 的前提下，所有非负分布里指数的熵最大（“固定均值”不能省——不加约束时熵没有上界）。用 Gibbs 不等式：对 \([0,\infty)\) 上任意均值为 \(\mu\) 的密度 \(p\)，取同均值的指数 \(q(x)=\tfrac1\mu e^{-x/\mu}\)，其 KL 散度非负（\(\int p\ln\frac{p}{q}\ge 0\)），展开：
（这是”猜答案再验证”：\(q\) 不是凭空取的——它正是 2.6 变分推出的候选；换个角度，它也是唯一能让下面交叉熵项塌成常数的 \(q\)，见《细节》。）
\[ \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} \]
移项即 \(H(p) \le 1 + \ln\mu = H(q)\)，等号仅当 \(p = q\)。固定均值下，指数是唯一的最大熵解——不是图里凑巧最高。
两条路的分工（贯穿全文的关键）：求最大熵分布分两步——2.6 的变分/拉格朗日负责”找出”候选 \(p\propto e^{-\lambda f}\)（但变分只给驻点，不自动保证是全局最大）；这里的 Gibbs 不等式负责”证明”它是唯一的全局最大——对所有同约束分布 \(H(p)\le H(q)\)，且等号仅当 \(p=q\)，所以任何别的分布都严格更小，指数是唯一达到最大的那个（代价是需要候选作输入）。一句话：变分找出它，Gibbs 盖章。
下面的《细节》解释这里用到的 KL 散度、以及为什么取指数作 \(q\)。只想要结论的话，可以直接跳到下一节《无记忆性》。
3.3.1 细节：KL 散度与两个选择
KL 散度是什么。 两个密度 \(p\)、\(q\) 之间的 KL 散度 \(D(p\|q) = \int p\ln\frac{p}{q}\,dx\)，读成”用 \(q\) 描述真实分布 \(p\) 时平均多付的对数代价”。这个”代价”有实在背景：香农编码给概率 \(p(x)\) 的符号分配 \(-\log_2 p(x)\) 比特，按真实 \(p\) 编码平均码长最短（\(=H(p)\)）；若误以为分布是 \(q\)、按 \(q\) 造码，每个符号平均多花的比特数恰好就是 \(D(p\|q)\)——用错模型永远不会更省，所以它 \(\ge 0\)。机器学习里的交叉熵损失正是这个量：训练分类器 = 让模型分布 \(q\) 逼近真实标签 \(p\)，即最小化 \(D(p\|q)\)。它永远非负，等号仅当 \(p=q\)：
自行观看（拓展）。 But what is cross-entropy? | Compression is Intelligence Part 2（3Blue1Brown）把交叉熵、编码代价与机器学习损失连成一条线；读完上面的编码解释再看最合适。
\[-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\]
那步 Jensen 不等式说的是：对凹函数（concave，∩ 形、开口朝下，如 \(\ln\)），“先平均、再取函数值” \(\ge\) “先取函数值、再平均”，即 \(\varphi(E[X]) \ge E[\varphi(X)]\)。几何上，凹函数的任意弦都在曲线下方，所以曲线在”平均点”处的高度高于两端函数值的平均。（注意别被汉字坑：字形”凹”像 ∪，但凹/concave 函数是 ∩；“凸”/convex 反而是 ∪ 形，如 \(x^2\)。盯例子别盯字。）\(\ln\) 处处凹（\(\ln'' x = -1/x^2 < 0\)），故 \(\ln E[Y] \ge E[\ln Y]\)；这里取 \(Y = q/p\)、期望按 \(p\) 算，把 \(\ln\) 提到期望外只会让值变大，得到 \(\le \ln E_p[q/p] = \ln\!\int q = \ln 1 = 0\)。于是 \(-D \le 0\)，即 \(D \ge 0\)。
为什么用 KL。 它恒 \(\ge 0\)，是把”\(=\)“撬成”\(\le\)“的现成工具。上面第一行拆开就得 \(H(p) \le -\int p\ln q\)，右边是交叉熵——熵不超过用任意 \(q\) 算的交叉熵。剩下的只是把交叉熵算出来。
为什么取指数。 一般的 \(q\) 交叉熵算不动。指数的 \(\ln q = -\ln\mu - x/\mu\) 对 \(x\) 线性，于是交叉熵里只出现 \(\int x\,p = p\) 的均值 \(=\mu\)（已固定）——它不再依赖 \(p\) 的其他细节，对每个均值 \(\mu\) 的 \(p\) 都塌成同一个常数 \(\ln\mu+1 = H(q)\)。起作用的正是 \(f(x)=x\) 的线性（2.6 通解里的 \(f\)）。
通用配方。 换约束换 \(q\)，同一个证明照走：约束 \(\int f\,p = c\) → 取 \(q \propto e^{-\lambda f}\)（\(\ln q\) 对 \(f\) 线性）→ 交叉熵只看 \(\int f\,p = c\) → 塌成常数。\(f=x\to\) 指数，\(f=x^2\to\) 正态，\(f=E\to\) 玻尔兹曼。这也说明”变分找出、Gibbs 盖章”这套分工对每个分布都成立，不止指数。
3.4 性质：无记忆性
指数分布还有一个独立于最大熵的特殊性质：唯一满足无记忆性的连续分布。
\[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")
数学本质：\(P(X > t) = e^{-\lambda t}\)——指数函数的”每增加一点、概率按固定比例下降”直接蕴含无记忆性。
3.5 物理实例：等温大气与放射性衰变
等温大气公式。一个空气分子在高度 \(h\) 处的重力势能是 \(E = mgh\)——能量随高度线性。这正是”线性约束 → 指数”的实物版，Boltzmann 因子直接给出：
\[P(h) \propto e^{-mgh/k_BT}\]
标高 \(H_{\text{scale}} = k_BT/mg \approx 8.5\) km（地球大气，\(T \approx 288\) K；大气物理里通常直接记作 \(H\)，本文加下标以免与熵的 \(H\) 撞名）：每升高 8.5 km，气压掉到 \(1/e\)。这也是《玻尔兹曼分布》那章的预告——同一个公式，那里换成能级。
放射性衰变。原子核不老化：一个已经存在了一万年的 \(^{14}\)C 核，接下来一秒衰变的概率和刚生成的核完全一样——这就是上一小节的无记忆性，所以寿命只能服从指数分布。半衰期 \(t_{1/2} = \ln 2/\lambda\) 等间距地把存活比例砍半：\(1/2, 1/4, 1/8, \dots\)。化学里的一级反应动力学（\(-d[A]/dt = k[A]\)）和荧光寿命衰减是同一个数学。
## ── 物理实例：等温大气 + 放射性衰变 ──────────────────────────
par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))
# 左图：等温大气公式 P(h) = exp(-h/H)，H = kT/mg ≈ 8.5 km
h_seq <- seq(0, 35, length.out = 400)
H_scale <- 8.5 # km，地球大气标高
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)
# 右图：放射性衰变，存活曲线 + 等间距半衰期
lam <- 1 # 衰变常数
t_half <- log(2) / lam
lifetimes <- rexp(N, rate = lam) # 模拟 N 个原子核的寿命
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 决定形状的是 \(f(x)\) 的形式
这两个例子给出指数分布，根源都在通解 \(p(x) \propto e^{-\lambda f(x)}\) 里 \(f(x)\) 的形式——\(\ln p\) 对 \(x\) 是几次函数，就得到哪个分布：
\[f(x) = x\ (\text{线性}) \to \text{指数}, \qquad f(x) = x^2 \to \text{半正态}, \qquad \dots\]
“只约束均值”写成数学就是 \(\int x\,p\,dx = \mu\)，即 \(f(x) = x\) 线性——这一步逼出指数。换一个被约束的量，形状就变：在同样的 \([0,\infty)\) 上若约束二阶矩（\(f(x) = x^2\)），得到的是半正态 \(p \propto e^{-\lambda x^2}\)。可见拍板的是被约束量对 \(x\) 的函数形式，而不是约束或参数的个数——单看”只有一个尺度参数”并不足以定出指数（Rayleigh、半正态也都是 \([0,\infty)\) 上的单参数分布）。
两个例子的线性各有来源：
- 等温大气：势能 \(E = mgh\) 对高度线性 → 指数。若势能 \(\propto h^2\)，同一个温度参数给出的会是高斯。
- 放射性衰变：对应的线性条件是恒定风险率（单位时间衰变概率不随年龄变），也就是无记忆性——\(\ln(\text{存活率}) = -\lambda t\) 对 \(t\) 线性。
一句话：指数 \(\Leftrightarrow \ln p\) 对 \(x\) 线性 \(\Leftrightarrow\) 恒定风险率 \(\Leftrightarrow\) 无记忆。指数、正态、玻尔兹曼的分别，就是 \(\ln p\) 各为 \(x\) 的一次、二次、能量函数——回到 2.6 节的通解。
4 幂律分布 = 最大熵 + 对数均值约束
幂律分布 \(p(x) \propto x^{-\lambda}\) 在自然界随处可见——地震能量、城市人口、词频、财富、网站链接数、物种多度……但课堂上很少讲它为什么出现，更少讲它和指数、正态的内在关系。答案出奇地简洁：同一个最大熵框架，只要换一个约束函数，就吐出幂律。前面指数约束的是算术均值 \(\langle x\rangle\)、正态约束的是方差，这里约束的是对数均值 \(\langle\ln x\rangle\)，等价于固定几何均值 \(e^{E[\ln X]}\)。
4.1 推导
约束：\(x \ge x_{\min} > 0\)，且对数的均值固定：
\[E[\ln X] = m \qquad (m > \ln x_{\min})\]
取约束函数 \(f(x) = \ln x\)，代入通解：
\[p(x) = \frac{1}{Z}\, e^{-\lambda \ln x} = \frac{1}{Z}\, x^{-\lambda}\]
这就是幂律。\(\ln p = -\lambda\ln x - \ln Z\) 对 \(\ln x\) 线性——所以幂律在双对数坐标上是一条直线，斜率 \(-\lambda\)。这条 log-log 直线正是”\(\ln p\) 对被约束量 \(\ln x\) 线性”的字面画面，和指数在 lin-log 上是直线完全平行。
归一化需要下限 \(x_{\min}\)。 幂律和指数一样，都不是全实轴上的自由分布：指数的自然支撑是 \([0,\infty)\)，幂律的自然支撑是 \([x_{\min},\infty)\)。差别在于，指数在 0 附近没有收敛问题，而 \(x^{-\lambda}\) 若一路延伸到 0 会发散，所以幂律必须有严格正的下限 \(x_{\min}>0\)。此时
\[Z = \int_{x_{\min}}^{\infty} x^{-\lambda}\,dx = \frac{x_{\min}^{\,1-\lambda}}{\lambda - 1} \quad(\text{需 }\lambda > 1\text{ 才收敛})\]
于是归一化的幂律（即 Pareto 分布）为
\[p(x) = (\lambda-1)\, x_{\min}^{\,\lambda-1}\, x^{-\lambda}, \qquad x \ge x_{\min}\]
记号约定：本文 \(\lambda\) 一律是密度指数（\(p\propto x^{-\lambda}\)，也就是最大熵通解里的那个乘子）。文献里常说的尾指数指生存函数 \(P(X>x)\propto x^{-(\lambda-1)}\) 的 \(\lambda-1\)——两者差 1。后文财富例子里的”生存指数 \(\alpha\)“就是 \(\lambda-1\)，别混。
\(\lambda > 1\) 是幂律”活得下去”的门槛（否则积不出有限概率），和指数要 \(x\ge0\)、正态无需截断，是同一类”支撑/收敛”前提。约束 \(m=E[\ln X]\) 进一步决定 \(\lambda\)：
\[ E[\ln X] = \ln x_{\min} + \frac{1}{\lambda - 1} \quad\Longrightarrow\quad \lambda = 1 + \frac{1}{m - \ln x_{\min}} \]
这一步很有解释力：如果 \(m\) 只比 \(\ln x_{\min}\) 大一点，说明尺度几乎贴着下限，\(\lambda\) 很大，尾巴很薄；如果 \(m\) 很大，说明几何均值远离下限，\(\lambda\) 逼近 1，尾巴极厚。
微分熵同样可以直接算：
\[ \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} \]
注意这里的 \(\lambda\) 是密度的幂指数：\(p(x)\propto x^{-\lambda}\)。因此 \(E[X]\) 有限需要 \(\lambda>2\)，方差有限需要 \(\lambda>3\)。\(\lambda>1\) 只保证”总概率能归一化”，不保证均值、方差都存在。
4.2 幂指数的三个门槛：1、2、3
幂律最容易误解的地方，是把”能归一化”、“均值有限”、“方差有限”混成一件事。它们其实是三个不同门槛：
| \(\lambda\) 区间 | 总概率 | 均值 \(E[X]\) | 方差 \(\mathrm{Var}(X)\) | 直觉 |
|---|---|---|---|---|
| \(\lambda \le 1\) | 不有限 | — | — | 不是合法概率分布 |
| \(1 < \lambda \le 2\) | 有限 | 无限 | 无限 | 极端大值强到连理论平均数都不存在 |
| \(2 < \lambda \le 3\) | 有限 | 有限 | 无限 | 有典型平均水平，但波动被极端大值支配 |
| \(\lambda > 3\) | 有限 | 有限 | 有限 | 仍是重尾，但平均和波动都可用常规矩描述 |
为什么是这些门槛？因为第 \(r\) 阶矩需要
\[E[X^r] = \int_{x_{\min}}^\infty x^r p(x)\,dx \propto \int_{x_{\min}}^\infty x^{r-\lambda}\,dx\]
而 \(\int^\infty x^a dx\) 只有在 \(a<-1\) 时才收敛，所以
\[E[X^r] < \infty \quad\Longleftrightarrow\quad r-\lambda < -1 \quad\Longleftrightarrow\quad \lambda > r+1\]
代 \(r=0\) 得 \(\lambda>1\)（总概率），代 \(r=1\) 得 \(\lambda>2\)（均值），代 \(r=2\) 得 \(\lambda>3\)（二阶矩/方差）。
## ── 幂指数门槛：归一化、均值、方差 ───────────────────────────
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)
区间 \(2<\lambda\le3\) 经常被单独讨论，因为它有一个很清楚的统计含义：平均规模存在，但波动没有有限的理论尺度。这时可以谈平均城市规模、平均财富、平均连接数；但方差不存在，样本波动会长期受少数极端大值支配。数据量增加时，新的超大观测仍可能显著改写方差估计。
另外两个区间也各有含义。\(1<\lambda\le2\) 更重尾：总概率已经有限，但理论均值不存在，样本平均会非常不稳定。\(\lambda>3\) 则温和得多：均值和方差都有限，虽然尾部仍比指数厚。实际数据还常有有限系统大小、测量下限、上限截断、机制混合，所以经验估计的 \(\lambda\) 可以落在不同区间；关键不是记某个固定范围，而是读懂 \(\lambda\) 所在区间对应的矩是否存在。
4.3 数值验证：固定 \(E[\ln X]\) 时，幂律的熵最高
取 \(x_{\min}=1\)，固定 \(E[\ln X]=0.5\)。最大熵解应为 \(\lambda = 1 + 1/0.5 = 3\) 的 Pareto 型幂律。下面把几个支撑同为 \([1,\infty)\)、且数值上调到同一个 \(E[\ln X]\) 的候选分布放在一起比较。
## ── 固定 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)
上面的每条候选分布都满足同一个约束；数值核验：
| 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 |
Power law 的红色条最长——在支撑 \([x_{\min},\infty)\) 和 \(E[\ln X]\) 固定时，幂律分布确实熵最高。
Uniform 看起来”平”，但它在 \(x>b\) 处把概率硬截成 0，偷偷加入了”绝不会超过 \(b\)“的额外信息；shifted exponential / gamma / Weibull 则把尾部压得比幂律薄，也是在约束之外多说了话。幂律只固定几何尺度，不固定算术均值，所以它允许极端大值出现得比指数多得多。
4.4 严格证明：固定对数均值下，幂律熵最大
这一节是可跳过的严格证明（与指数分布那节的 Gibbs 论证同构）。只想要结论——“固定 \(E[\ln X]\) 时幂律熵最大”——的话，可以直接跳到下一节《scale-free》。
用与指数分布完全同构的 Gibbs 不等式。对 \([x_{\min},\infty)\) 上任意满足 \(E_p[\ln X]=m\) 的密度 \(p\)，取同约束的幂律候选
\[q(x) = (\lambda-1)x_{\min}^{\lambda-1}x^{-\lambda}, \qquad \lambda = 1 + \frac{1}{m-\ln x_{\min}}\]
KL 散度非负：
\[ \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} \]
最后一行正是 \(H(q)\)，所以 \(H(p)\le H(q)\)，等号仅当 \(p=q\)。这就是”变分找出它，Gibbs 盖章”的同一套路：因为 \(\ln q\) 对 \(\ln x\) 线性，交叉熵只依赖固定的 \(E[\ln X]\)，不再依赖 \(p\) 的其他细节。
4.5 幂律 = 对数尺度上的指数
这一节把”幂律”彻底还原成一句话：它就是指数分布，只不过住在对数尺度上。先从几何看见它，再一步换元坐实。
几何入口：拉直它的坐标系，就读出它的约束。 通解取对数，\(p\propto e^{-\lambda f(x)}\) 两边取 \(\ln\)：
\[\ln p(x) = -\lambda\, f(x) + \text{常数}\]
\(\ln p\) 是 \(f(x)\) 的一次函数。把纵轴取 \(\ln p\)、横轴取 \(f(x)\)，任何最大熵分布都是一条斜率 \(-\lambda\) 的直线。反过来，“哪种坐标系把分布拉直” = “它约束哪个 \(f\)” = “它是哪种最大熵分布”——log-log、lin-log 不是随手挑的画法，是”把 \(\ln p\) 对被约束量画出来”。
## 每个最大熵分布只在"横轴 = 它的 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: 横轴 x → 指数直
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: 横轴 ln x → 幂律直
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: 横轴 x^2 → 高斯直
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)
三张图同一批曲线，只换横轴：每个分布只在”横轴 = 它自己的 \(f\)“那张里才笔直。读成对照表：
| 分布 | 把它拉直的坐标（纵 \(\ln p\)，横 \(\cdots\)） | 斜率 | 约束 \(f(x)\) |
|---|---|---|---|
| 指数 | \(x\)（lin-log） | \(-\lambda\) | \(x\) |
| 幂律 | \(\ln x\)（log-log） | \(-\lambda\) | \(\ln x\) |
| 正态 | \(x^2\) | \(-\tfrac{1}{2\sigma^2}\) | \(x^2\) |
点破：幂律的 log-log 直线，就是它在对数尺度上是指数。 幂律在 log-log 上直，意思是 \(\ln p\) 对 \(\ln x\) 线性——而”\(\ln p\) 对某个变量线性”正是指数分布在那个变量上的招牌（指数在 lin-log 上直）。所以幂律的 log-log 直线，等价于”把 \(\ln x\) 当新变量后，它是一条 lin-log 直线 = 一个指数分布”。换元把它坐实：令 \(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})\) 恰好服从指数分布，速率 \(\lambda-1\)：
\[\text{幂律}(x)\;\Longleftrightarrow\;\text{指数}(\ln x)\]
“\(\ln x\) 之于幂律 = \(x\) 之于指数”。这也顺手解释了推导里的约束公式 \(E[\ln X]=\ln x_{\min}+\tfrac{1}{\lambda-1}\)——右边的 \(\tfrac{1}{\lambda-1}\) 正是速率 \(\lambda-1\) 的指数分布的均值。约束几何均值 \(\langle\ln x\rangle\)，就是约束”对数空间里那个指数”的均值，和指数约束 \(\langle x\rangle\) 一模一样。所以幂律不是”另一种”指数，它就是指数，只是换了尺度。
由此，无记忆性的对仗也现形了：
| 指数分布 | 幂律分布 | |
|---|---|---|
| 无记忆形式 | \(P(X>s{+}t\mid X>s)=P(X>t)\) | \(P(X>ct\mid X>t)=P\!\left(\tfrac{X}{x_{\min}}>c\right)\) |
| 不变的操作 | 加一个量 \(t\)（平移） | 乘一个倍数 \(c\)（缩放） |
| 特征尺度 | 有，\(1/\lambda\) | 无（scale-free） |
指数是加性无记忆（线性尺度上平移不变），幂律是乘性无记忆（对数尺度上平移不变 = 缩放不变）。取一次对数，乘性就变加性、幂律就变指数——同一个”无记忆”，只是活在不同尺度上。下一节把这个”乘性无记忆”画出来，并证明幂律是唯一的无标度分布。
4.6 Scale-free：把”乘性无记忆”画出来，并证明它唯一
上一节已得幂律是乘性无记忆——\(P(X>ct\mid X>t)=c^{-(\lambda-1)}\)，只看倍数 \(c\)、不看起点 \(t\)。这里做两件上节没做的事：把它画出来，并证明幂律是唯一这样的分布。
## ── 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)
左图的水平线就是无标度：对幂律，“再翻倍”的概率不随起点 \(t\) 改变；右图指数则越往后再翻倍越难（\(P(X>ct\mid X>t)=e^{-(c-1)t/\mu}\) 依赖 \(t\)）——这正是指数有特征尺度 \(\mu\)、幂律没有的根源。
唯一性：若一个足够规则的正函数对任意缩放都满足 \(f(cx)=a(c)f(x)\)，它只能是 \(Kx^{-\lambda}\)。所以”精确无标度”这一条本身就把分布钉死成幂律——scale-free 不是幂律的一条性质，而是它的另一个定义。（经验数据只在尾部或有限区间近似满足；截断幂律、对数正态在有限 log-log 图上也近似直，所以直线是线索、不是证明，严格判定见《怎样才算”真的”幂律》。）
4.7 无标度 = 分形 = 临界：同一枚硬币的三面
刚证的”精确无标度只能是幂律”，一脚就踏进了分形几何——因为”无标度”正是分形的定义。这一节把幂律、分形、临界动力学接成同一件事，也是全章通向复杂系统的出口。
自相似 ⟺ 尺度不变 ⟺ 幂律。 分形的定义性质是自相似：放大任意倍数看起来一样。“缩放后形状不变”就是 \(f(cx)=a(c)f(x)\)——正是上面唯一性定理的前提，而结论是它只能是 \(x^{-\lambda}\)。所以这三句话是同一个命题的三种说法：分形自相似、尺度不变、幂律。
分形维数就是一个幂指数。 怎么测一个分形？用边长 \(\varepsilon\) 的盒子盖住它、数盒子数 \(N(\varepsilon)\)，分形满足 \(N(\varepsilon)\propto\varepsilon^{-D}\)——一条幂律，指数 \(D\) 就是分形维数。测分形 = 在 log-log 上量斜率（还是上一节那套）：
## Koch 曲线 + 盒计数：分形维数 D 就是一条 log-log 幂律的斜率
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)
Koch 曲线盒计数的斜率给出 \(D=\ln4/\ln3\approx1.26\)——非整数维，正是幂指数的化身。Cantor 集 \(\ln2/\ln3\approx0.63\)、Sierpinski \(\ln3/\ln2\approx1.58\)、海岸线 \(L(\varepsilon)\propto\varepsilon^{1-D}\)（量得越细越长）同理。所以幂律的指数 \(\lambda\) 和分形维数 \(D\) 是同一类东西——“某个量随尺度怎么标度”的指数。
动力学来路：混沌、SOC、临界。 幂律/分形在自然界满地都是，因为好几类动力学天然收敛到尺度不变：确定性混沌落在分形维的奇怪吸引子上（Lorenz、Hénon）；自组织临界（Bak–Tang–Wiesenfeld 沙堆）让耗散系统自动停在临界点、雪崩大小服从幂律，不用调参——地震（Gutenberg–Richter，就是前面表里那条）、森林火灾、神经雪崩都属此类；相变临界点上关联长度发散、没有特征尺度、关联按幂律衰减，重整化群（RG）把系统看成缩放变换的不动点，这是”尺度不变 \(\Rightarrow\) 幂律”最深的理论。
焊接（最漂亮的一处）。 本章给过幂律一条最大熵来路：在尺度不变测度 \(dx/x\) 下”什么都不知道”（见后面延伸章）。而 \(dx/x\) 正是唯一在 \(x\to cx\) 下不变的测度——它就是 RG／分形语言里的”尺度不变”。于是
\[\underbrace{\text{最大熵：对数尺度上最无知}}_{\text{推断语言}}\;=\;\underbrace{\text{RG：缩放变换的不动点}}_{\text{动力学语言}}\]
两种语言，同一个”没有特征尺度”。上面那条唯一性定理是它的静态影子，RG 是它的动力学版本。
给生物的一条线索。 分形/无标度在生物学里是结构性的：血管、支气管、神经树都有分枝网络；West–Brown–Enquist 曾从空间填充的分枝输送网络推导出 \(3/4\) 代谢标度，神经雪崩则促成了”临界脑”假说。(West et al. 1997; Beggs and Plenz 2003) 这些是有力的机制线索，而不是“所有生命系统都处在临界点”的定理。
4.8 合体：同时约束 \(\langle x\rangle\) 和 \(\langle\ln x\rangle\) → Gamma
指数只固定算术均值 \(\langle x\rangle\)，幂律只固定几何均值 \(\langle\ln x\rangle\)。两个同时固定呢？两个约束配两个乘子 \(\lambda_1\)（配 \(x\)）、\(\lambda_2\)（配 \(\ln x\)），代入通解：
\[p(x)\propto e^{-\lambda_1 x - \lambda_2\ln x} = x^{-\lambda_2}\,e^{-\lambda_1 x}\]
幂律核 \(\times\) 指数尾。在 \((0,\infty)\) 上这正是 Gamma 分布（形状 \(k=1-\lambda_2\)、速率 \(\lambda_1\)）；在 \([x_{\min},\infty)\) 上就是文献里的”带指数截断的幂律”。两个旋钮，两个极限：
- \(\lambda_1\to 0\)（松开算术均值）→ 指数尾消失 → 纯幂律 \(x^{-\lambda_2}\)；
- \(\lambda_2\to 0\)（松开几何均值）→ 幂律核消失 → 纯指数 \(e^{-\lambda_1 x}\)。
所以指数和幂律不是两个物种，是同一个族的两端：\(\lambda_2\) 管小-中尺度的幂律核（log-log 斜率 \(-\lambda_2\)），\(\lambda_1\) 管大尺度的指数截断（在 \(x\sim 1/\lambda_1\) 处把尾巴掰下来）。下面两节的细胞谱系次临界情形、财富的 Kesten 稳态，尾部都恰好是这个”核 + 截断”形——它们不是巧合地”像” Gamma，而就是两个约束下的最大熵解。这也是”双约束”的第一次登场：下一章正态把约束换成 \(x\) 和 \(x^2\)，同样两个乘子。
## Gamma 桥：p ∝ x^{-λ2} e^{-λ1 x}。两个旋钮：λ1=指数尾截断，λ2=幂律核斜率
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))
# 左：固定核 λ2=2，变截断 λ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)
# 右：固定截断 λ1=0.01，变核 λ2；λ2=0 = 纯指数
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)
左图：固定核 \(\lambda_2\)、增大 \(\lambda_1\)，指数截断的”膝盖”从右往左移，\(\lambda_1=0\) 时退回一条笔直的纯幂律。右图：固定截断 \(\lambda_1\)、增大核 \(\lambda_2\)，log-log 上的核越来越陡，\(\lambda_2=0\) 时退回纯指数。现实里”有限尺寸截断的幂律”几乎都活在这张图内部。
4.9 现实里的幂指数落在哪
前面都是框架，这里落地。一些被反复研究的幂律，按尾指数 \(\alpha\)（生存函数 \(P(X>x)\propto x^{-\alpha}\)，即密度指数 \(\lambda-1\)）列出——落在哪个矩区间（区间含义见前面《幂指数的三个门槛》），直接决定”平均、方差还有没有意义”：
| 现象 | 量 | 尾指数 \(\alpha\ (\approx)\) | 密度 \(\lambda=\alpha{+}1\) | 矩状况 |
|---|---|---|---|---|
| Zipf 词频 / 城市人口 | 频次、人口 | \(\sim 1\) | \(\sim 2\) | 均值临界、方差无限 |
| 财富 / 收入上尾（Pareto） | 财富 | \(1\text{–}2\) | \(2\text{–}3\) | 常在”均值有限、方差无限” |
| 网络度分布（无标度网） | 节点连边数 | \(2\text{–}3\) | \(3\text{–}4\) | 均值有限，方差常无限 |
| 论文引用 / 网站入链 | 被引、入链数 | \(2\text{–}3\) | \(3\text{–}4\) | 同上 |
| 地震（Gutenberg–Richter） | 释放能量 | \(\sim 2/3\) | \(\sim 5/3\) | 能量均值都不收敛，极端事件主导 |
| 基因家族大小（基因组演化） | 同源家族成员数 | \(\sim 2\) | \(\sim 3\) | 基因重复 + 分化（优先增长） |
一句话提醒：这些指数都是经验估计、区间近似。严格统计检验表明，很多被称作”幂律”的数据在整段上并不是干净幂律，而是截断幂律、对数正态或混合。(Clauset et al. 2009) 所以看到 log-log 近似直线，是线索，值得追问机制；不是终点。这个”线索不是证明”到底有多大坑，下一节展开。
4.10 怎样才算”真的”幂律：拟合、检验与陷阱
如果你要对真实数据（基因表达尾、网络度数……）判断是不是幂律，这一节是必读；只想懂原理的话可以跳到《现实案例》。核心参考是 Clauset、Shalizi、Newman 的方法论综述。(Clauset et al. 2009)
两个常见做法都是错的。 早期大量”幂律”论文只靠：(1) 把数据画在 log-log 上、肉眼看直线；(2) 对 log-log 直方图做最小二乘线性回归、报一个高 \(R^2\)。两条都不成立：
- log-log 回归的斜率估计有偏——分箱假象、尾部箱子噪声大、最小二乘的误差假设在对数空间根本不成立；
- 高 \(R^2\) 证明不了幂律——对数正态、指数、拉伸指数在有限区间的 log-log 上都能”看着挺直”。眼球更不是检验。
正确的四步配方（现在的金标准）：
- 先定下限 \(x_{\min}\)：幂律通常只在尾部成立。对每个候选 \(x_{\min}\)，用 Kolmogorov–Smirnov 距离（数据 vs 拟合幂律）最小化来选它。
- 用最大似然估 \(\lambda\)，不用回归。连续情形就是 Hill 估计量 \(\hat\lambda = 1 + n\big/\sum_i \ln\frac{x_i}{x_{\min}}\)，相合无偏；回归有偏。
- 拟合优度检验：算 KS 统计量，再用自助重采样得 \(p\) 值。\(p\) 小 → 拒绝幂律假设。（\(p\) 大只是”没拒绝”，不等于”证明是幂律”。）
- 模型比较：用似然比检验把幂律和对数正态、指数、拉伸指数、带指数截断的幂律对打。真问题不是”是不是幂律”，而是”幂律比这些替代品更好吗”。
结论很扎心。 Clauset 等把 24 个前人称作幂律的真实数据集重估了一遍：干净通过的只有少数；很多和对数正态无法区分，或带指数截断的幂律拟合更好。后续 Broido 与 Clauset 重审约一千个网络，题目直接是《Scale-free networks are rare》。(Clauset et al. 2009; Broido and Clauset 2019)
是、还是不是？对照才有意思。 同样”看着重尾”，有的是幂律，有的根本不是——而判据回到过程（本章一直在讲的那条轴）：
| 现象 | 常被说成 | 更站得住的 | 生成过程 |
|---|---|---|---|
| 词频（Zipf）、地震能量 | 幂律 | 幂律（支持较强） | 临界 / 优先增长 |
| 财富、收入的主体 | 幂律 | 对数正态（仅极上尾近 Pareto） | 乘性增长，无托底 |
| 多数网络度分布 | 无标度幂律 | 常带截断或对数正态；纯无标度罕见 | —（Broido–Clauset 2019） |
| 身高、体重、测量误差 | —— | 高斯 | 加性平均 → CLT |
| 放射性寿命、无记忆等待 | —— | 指数 | 恒定风险率 |
四种过程，四种命运：加很多小量 → 高斯；无记忆等待 → 指数；乘很多小量、没有下界托底 → 对数正态；乘性 + reset/托底（Kesten）或临界/优先连接 → 真幂律。所以”是不是幂律”这个经验问题，本质是”背后是哪种过程”。重尾 \(\ne\) 幂律：对数正态同样重尾，却来自完全不同的机制，在有限数据上还常和幂律难分——这正是 Clauset 检验要拆开的核心，也是为什么”看到直线要追问机制、而不是宣布胜利”。
这不是说”幂律不存在”，而是：证据要用 MLE + 拟合优度 +
模型比较,并老实承认替代分布常常拟合得一样好、甚至更好。注意它点名的两个头号劲敌——对数正态和带指数截断的幂律——正是本章已经出现过的两位：前者来自乘性过程
+ CLT（《旁例：财富》），后者就是 Gamma
桥的内部。换句话说，“到底是幂律还是它俩”这个经验难题，本章的最大熵框架早已把三者摆在同一张谱系图上。实操工具：R
的 poweRlaw 包、Python 的 powerlaw
包都实现了上面这套配方。
上面那张过程表是”总”。接下来《现实案例》把其中两行做实——临界分枝（产出真幂律 \(T^{-3/2}\)）与乘性 + 托底的 Kesten 过程（主体对数正态、上尾 Pareto、常带截断）——看具体机制怎样产出幂律，也顺便看清那两个”劲敌”是怎么从同一批机制里冒出来的。
4.11 现实案例：细胞谱系的临界分枝怎样给出幂律尾部
这是上一节过程分类表里的第一行（临界 / 分枝）做实：幂律不只能由”只知道几何尺度”的最大熵推出，也能从一个具体的细胞分裂机制里长出来。本节结论一句话：恰好临界（分裂=终止）的分枝过程，累计产出服从 \(T^{-3/2}\) 幂律；稍微偏离临界，尾部就被指数截断（这就落进了 Gamma 桥的”幂律核 + 指数尾”）。这是 Galton–Watson 分枝过程的经典结果。(Harris 1963) 想看这个结论怎么从组合数精确算出来，往下读；只要记住这句话的话，可以跳到《旁例：财富》。
考虑一个刻意简化的细胞谱系模型：从一个创始细胞开始，每个细胞独立地以概率 \(q\) 分裂成两个子细胞，以概率 \(1-q\) 终止（可理解为死亡或终末分化）；两个子细胞再重复同一规则。
这里研究的不是某个时刻的克隆大小，而是直到谱系终止为止，累计出现过的总细胞数 \(T\)，创始细胞也算一个。若发生 \(n\) 次分裂，谱系恰有 \(n+1\) 次终止，因而 \(T=2n+1\)。满足这一条件的二叉谱系形状有第 \(n\) 个 Catalan 数
\[ C_n=\frac{1}{n+1}\binom{2n}{n} \]
种；每一种形状有同样的概率。因此
\[ P(T=2n+1)=C_nq^n(1-q)^{n+1}. \]
Catalan 数在这里并不神秘：它数的是“分裂产生的活细胞从未先于终止事件耗尽”的所有合法谱系树。这个封闭形式让我们能直接看见临界点发生了什么。
当 \(q=1/2\)，分裂与终止恰好平衡，是临界分枝。用 \(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}. \]
所以这个离散机制产生了幂指数 \(3/2\) 的渐近尾部。它落在前面的 \(1<\lambda\le2\) 区间：理想的、无限持续的临界模型中，\(E[T]\) 不存在。不是说一片真实组织会产生无限多细胞，而是说极少数极长谱系足以让理论平均值失去稳定尺度。
## ── 细胞谱系：临界分枝给出幂律，偏离临界即出现截断 ─────────
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)
左图中，\(q=1/2\) 的红线在 log-log 坐标下趋近直线；\(q=0.45\) 或 \(0.40\) 时，尾部向下弯。右图把 \(T^{-3/2}\) 这一幂律因子除掉：临界曲线趋向常数 \(\sqrt{2/\pi}\)，而两条次临界曲线继续掉到 0。这不是肉眼拟合，而是上面 Catalan 公式的精确计算。
更精确地，\(q<1/2\) 时
\[ P(T=t)\asymp t^{-3/2}\exp\left[-\frac{t}{2}\log\frac{1}{4q(1-q)}\right]. \]
也就是“幂律核 \(\times\) 指数截断”。\(q\) 离 \(1/2\) 越远，截断越早出现；若 \(q>1/2\)，则有正概率永不终止，\(T\) 不再是一个只取有限值的概率分布。真实组织还会受到空间、资源、发育时间窗口和调控反馈的限制，因此即使局部接近临界，观测到的尾部也会被进一步截断。
这个案例与前文的最大熵推导回答的不是同一个问题。前文的连续 Pareto 是：在 \([x_{\min},\infty)\) 上只固定 \(E[\ln X]\) 时，最少额外假设的分布；这里是：给定一个细胞分裂/终止的微观规则，累计谱系产出会怎样。前者给出推断原则，后者给出一种产生近似幂律的细胞机制。两者并不互相替代。
4.11.1 旁例：财富分布为什么常被谈幂律？
这是过程分类表的第二、三行做实（乘性增长、以及乘性 + 托底）。财富是个很有用的对照，恰好纠正一句常见错话：“乘法增长必然产生幂律”并不对。 若只有比例增长
\[ W_{t+1}=A_tW_t, \]
那么 \(\ln W_t=\ln W_0+\sum_{s<t}\ln A_s\)；许多小的独立增长率相加时，中心极限定理首先指向的是对数正态，而不是幂律。
产生 Pareto 尾的一类经典机制是在比例增长之外持续加入收入、最低保障或重新进入：
\[ W_{t+1}=A_tW_t+B_t,\qquad B_t>0. \]
在适当条件下（例如 \(E[\ln A]<0\) 使总体过程有稳态，但仍有一部分时期 \(A>1\)），这个 Kesten 型递推的上尾满足
\[ P(W>w)\sim w^{-\alpha}, \qquad E[A^\alpha]=1, \]
其中 \(\alpha\) 是尾/生存指数（对应前面记号约定的 \(\lambda-1\)，即密度指数 \(\lambda=\alpha+1\)）。直觉上，\(B_t\) 不让财富掉到零，偶尔连续的正增长又能把少数轨迹推得很高；两者共同形成”主体有尺度、上尾近似无尺度”的稳态。(Kesten 1973; Gabaix 2009)
所以财富数据里常见的不是一条从底到顶的纯幂律，而可能是主体近似对数正态、较高阈值以上才接近 Pareto 尾；税收、破产、有限寿命和财富上限都会让尾部弯下去。它与细胞谱系例子共同说明：幂律不是一个万能标签，而是由具体机制、观察尺度和截断条件共同决定的尾部描述。
5 正态分布 = 最大熵 + 均值 + 方差约束
5.1 推导
约束：\(E[X] = \mu\)，\(\text{Var}(X) = \sigma^2\)。需要两个约束函数 \(f_1(x) = x\) 和 \(f_2(x) = x^2\)，代入通解：
\[p(x) = \frac{1}{Z}\, e^{-\lambda_1 x - \lambda_2 x^2}\]
配方：指数里是 \(x\) 的二次式，对它配方（把一次项凑进平方），常数项并入 \(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{并入 }Z}\]
于是 \(p(x) \propto e^{-\lambda_2 (x - \mu)^2}\)，是中心在 \(\mu = -\tfrac{\lambda_1}{2\lambda_2}\)、方差 \(\tfrac{1}{2\lambda_2}\) 的高斯。对照约束得 \(\lambda_2 = \tfrac{1}{2\sigma^2}\)、\(\lambda_1 = -\tfrac{\mu}{\sigma^2}\)，归一化常数 \(Z = \sqrt{2\pi\sigma^2}\)：
\[p(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\, e^{-(x - \mu)^2 / (2\sigma^2)}\]
微分熵 \(H = -\int p\ln p\,dx\)（步骤与 3.1 同构，只是把”均值”换成”方差”）：
\[ \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{拆开} \\[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} \]
即 \(H[\text{Normal}(\mu,\sigma^2)] = \tfrac12\ln(2\pi e\,\sigma^2)\)。注意它只依赖 \(\sigma^2\)、与 \(\mu\) 无关——平移分布不改变展开程度。
自行观看（拓展）。 Why π is in the normal distribution (beyond integral tricks)（3Blue1Brown）专门解释上面归一化常数里的 \(\pi\) 从何而来；适合在接受高斯形式之后观看。
5.2 数值验证：候选分布中，正态的熵最高
在 \(\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
}
# 均值 = 0，方差 = 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)
Normal 的红色条最长——在均值和方差都固定时，正态分布的熵确实最大。
Laplace 和 Logistic 尽管形状接近，但尾部行为不同——要么太重（Laplace），要么太轻（Uniform），都额外注入了假设。
5.3 CLT：另一条路到同一个终点
中心极限定理说的是：无论原始分布是什么，大量独立样本的均值趋向正态。
自行观看。 But what is the Central Limit Theorem?（3Blue1Brown）以动画展示“不同形状为什么会在求和后趋同”；可先看视频再读下面的模拟图，也可读完后对照。
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 和最大熵是到达正态分布的两条独立道路：
- CLT 说：加法的极限是正态（数学事实）
- 最大熵说：在固定均值和方差下，最无偏见的分布是正态（逻辑原则）
它们的结论一致，但逻辑完全独立。CLT 不需要熵的概念，最大熵不需要求和的概念。
下面用傅里叶变换给出 CLT 的证明。只想要结论的话，可以直接跳到下一节《物理实例》。
5.3.1 细节：用傅里叶变换证明 CLT
思路。 上一段说过：独立随机变量相加，密度做卷积。\(n\) 个相加就是密度自己卷 \(n\) 次——直接算是噩梦。但傅里叶变换有个招牌性质：把卷积变成乘法。于是”卷 \(n\) 次”变成”乘 \(n\) 次方”，立刻好算。这就是整个证明的引擎。
自行观看（前置补充）。 But what is a convolution?（3Blue1Brown）解释卷积的几何图像。它不是理解前面 CLT 结论的前提，但会让下面的傅里叶证明更顺。
要证什么。 \(X_1,\dots,X_n\) 独立同分布，均值 \(\mu\)、方差 \(\sigma^2\)（有限）。标准化的和 \(Z_n = \dfrac{\sum_i X_i - n\mu}{\sigma\sqrt{n}}\) 依分布收敛到 \(N(0,1)\)。先中心化标准化 \(Y_i = (X_i-\mu)/\sigma\)（\(E[Y_i]=0\)、\(E[Y_i^2]=1\)），则 \(Z_n = \tfrac{1}{\sqrt n}\sum_i Y_i\)。
工具：密度的傅里叶变换 \(\hat f(t) = \int f(x)\,e^{itx}\,dx = E[e^{itX}]\)（概率里叫它特征函数，本质就是密度的傅里叶变换）。要用两条傅里叶的招牌性质，外加一个收敛桥：
- 卷积定理：\(\widehat{f * g} = \hat f\cdot\hat g\)（时域卷积 = 频域相乘）。所以独立和的变换 = 各自变换相乘——这一条把”和”变成”积”。
- 高斯自变换：\(N(0,1)\) 的密度的傅里叶变换是 \(e^{-t^2/2}\)（高斯变换后还是高斯）。
- Lévy 连续性定理（收敛桥，当黑箱）：\(\hat f\) 逐点收敛 \(\Rightarrow\) 分布收敛。
计算。 由卷积定理，\(\hat f_{Z_n}(t) = \Bigl[\hat f_Y\!\bigl(\tfrac{t}{\sqrt n}\bigr)\Bigr]^n\)。把 \(\hat f_Y\) 在 0 处 Taylor 展开（这一步用到均值 0、方差 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{展开 }e^{isY} \\[2pt] &= 1 - \tfrac{s^2}{2} + o(s^2) && E[Y]=0,\ E[Y^2]=1 \end{aligned} \]
代 \(s = t/\sqrt n\)，得 \(\hat f_Y(\tfrac{t}{\sqrt n}) = 1 - \tfrac{t^2}{2n} + o(\tfrac1n)\)，于是
\[\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}\]
（用 \((1+\tfrac an+o(\tfrac1n))^n \to e^a\)，\(a=-t^2/2\)。）
到这里只证明了变换逐点收敛 \(\hat f_{Z_n}(t)\to e^{-t^2/2}\)，还差最后一步把它变成分布收敛。\(e^{-t^2/2}\) 恰是 \(N(0,1)\) 的变换（性质 2，负责”认出”是哪个分布）；再由 Lévy 连续性定理（工具 3）——特征函数逐点收敛、且极限在 \(t=0\) 连续（\(e^0=1\) ✓），则分布收敛——得 \(Z_n \xrightarrow{d} N(0,1)\)。∎（\(t=0\) 连续这一条不是摆设：它排除概率质量逃向无穷的退化情形。）
为什么偏偏是高斯。 看展开 \(\hat f_Y(s) = 1 - \tfrac{s^2}{2} + \cdots\)：\(\tfrac{1}{\sqrt n}\) 这个缩放精确调准，让二阶项（方差）在 \(n\) 次幂后恰好活下来、收成 \(e^{-t^2/2}\)，而三阶及以上（偏度……）带更高次的 \(1/\sqrt n\)，被碾成零。极限里只有”均值 0、方差 1”幸存，其余全抹掉——所以不管原分布长什么样（只要方差有限），和都收敛到同一个高斯。这与最大熵正好对上：固定均值和方差、其余不管 \(\to\) 正态。两个彩蛋：(1) \(e^{-t^2/2}\) 本身就是高斯，指数上的二次与正态 \(\ln p\) 的二次是一回事；(2) 高斯是卷积的不动点（两个高斯卷积还是高斯），所以极限稳定停在高斯、不再变化。
5.4 物理实例：化学键振动与布朗扩散
化学键振动。把化学键近似成弹簧，偏离平衡位置 \(x\) 的势能是二次的：\(E = \frac{1}{2}\kappa x^2\)。泡在温度 \(T\) 的热浴里，Boltzmann 因子作用到二次能量上：
\[p(x) \propto e^{-\kappa x^2 / 2k_BT} = \text{Normal}\!\left(0,\ \sigma^2 = k_BT/\kappa\right)\]
“二次约束 → 正态”的分子级实例：原子在平衡键长附近做高斯分布的热涨落，温度越高分布越宽（\(\sigma \propto \sqrt{T}\)），振动光谱里直接可测。气体分子的速度分量是同一个逻辑（动能 \(\frac{1}{2}mv_x^2\) 也是二次的），见《玻尔兹曼分布》那章。
这里为什么得到高斯，而不是”玻尔兹曼分布”？ 因为玻尔兹曼因子是机制，具名分布由 \(E(x)\) 的次数决定。\(e^{-E/k_BT}\) 只说”概率随能量指数衰减”，没说坐标 \(x\) 服从什么——那取决于能量怎么依赖 \(x\)：
- 势能对坐标一次（等温大气的 \(E = mgh\)）→ \(e^{-mgh/k_BT}\) 是 \(h\) 的指数分布；
- 势能对坐标二次（这里的 \(E = \frac12\kappa x^2\)，或动能 \(\frac12 mv_x^2\)）→ \(e^{-\kappa x^2/2k_BT}\) 是 \(x\) 的高斯分布。
同一个因子，换一个 \(E(x)\)，就换一个具名分布。所以指数与高斯不是玻尔兹曼分布的”竞争者”，而是它在不同势能形式下的两种面貌——下一章会把这条主线连同 Maxwell 速率分布一起走完。
布朗扩散。墨水滴进水里，每个色素粒子的位移是无数次分子碰撞冲量的累加——这正是上一小节 CLT 的物理化身（Einstein 1905）。\(t\) 时刻的位置分布是方差 \(\propto t\) 的高斯，剖面随 \(\sqrt{t}\) 展宽。
## ── 物理实例：化学键振动 + 布朗扩散 ──────────────────────────
par(mfrow = c(1, 2), mar = c(4.6, 5.0, 3.5, 1))
# 左图：谐振子势 E = kappa x^2 / 2 泡在热浴里 → 位置分布是正态
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)
# 右图：布朗扩散——随机行走的位置分布是变宽的高斯
n_part <- 20000
t_snap <- c(25, 100, 400)
cols_d <- c("steelblue", "purple", "tomato")
# 每个粒子 = 许多次独立碰撞冲量的累加（每步 ~ 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 玻尔兹曼分布 = 最大熵 + 能量约束
6.1 Boltzmann 原理就是最大熵原理
热力学系统在温度 \(T\) 下：
\[\boxed{P(\text{state}) \propto e^{-E / k_B T}}\]
这个公式常被当作一个”物理定律”来记忆。但它的推导就是最大熵：
- 约束：总能量守恒，\(\langle E \rangle\) 固定
- 最大熵解：\(p \propto e^{-\lambda E}\)
先说清这一章与前几章的层级关系。 开篇那张表里，前三行的约束是 \(x\) 的具体函数（\(x\)、\(\ln x\)、\(x^2\)），而这一行写的是 \(f(x) = E(x)\)——它是个总纲，不是第四个并列的具名分布。\(E(x)\) 一旦落地，掉出来的就是前面已经见过的那几个：
| 能量怎么依赖变量 | \(e^{-E/k_BT}\) 成为 | 实例 |
|---|---|---|
| \(E\) 对变量一次：\(E = mgh\)、等间距能级 | 该变量的指数分布 | 等温大气、能级占据 |
| \(E\) 对变量二次：\(E = \frac12 mv_x^2\)、\(\frac12\kappa x^2\) | 该变量的高斯分布 | 速度分量、化学键振动 |
| 二次能量 \(+\) 非均匀态密度 \(g(v)\propto v^2\) | Maxwell 速率分布 | 分子速率 |
所以”指数、高斯、玻尔兹曼”不是三个互相竞争的分布：玻尔兹曼因子是共用的机制，指数和高斯是它在一次 / 二次能量下的两张脸，Maxwell 速率则是再乘上态密度的结果。本章按这个顺序逐个走一遍。
顺着这一层，得说清一件容易被高估的事：玻尔兹曼形式本身没有排除任何分布。 给定任意一个处处为正的密度 \(p(x)\)，只要定义
\[E(x) \;\equiv\; -k_BT\ln p(x)\]
立刻就有 \(p(x)\propto e^{-E(x)/k_BT}\)。也就是说，“这是玻尔兹曼分布”作为一句关于变量 \(x\) 的主张，是同义反复——任何正密度都写得成这个样子。
换句话说，\(e^{-E/k_BT}\) 这个形式本身不含信息量。玻尔兹曼因子的全部力量，在于物理独立地告诉你 \(E(x)\) 是什么——重力给 \(mgh\)、弹簧给 \(\frac12\kappa x^2\)、动能给 \(\frac12 mv^2\)。\(E\) 是从外部输入的，不是从数据里拟合出来的；一旦 \(E\) 被钉死，\(e^{-E/k_BT}\) 就变回一个只剩温度这一个参数的、可以被证伪的分布。
本章接下来就沿着”\(E\) 被物理钉死之后能榨出多少东西”往下走：配分函数如何把温度的一次求导变成能量、二次求导变成涨落；能级的简并如何在玻尔兹曼因子之外再乘一层；以及为什么每个二次自由度恰好分到 \(\frac12 k_BT\)。
\(\lambda\) 就是拉格朗日乘子。物理学家给它取了个名字：\(\beta = 1/(k_BT)\)。
也就是说：温度不是一个独立的物理量，而是最大熵推导中拉格朗日乘子的倒数。
- 高温（\(\beta\) 小）→ 约束松 → 分布更展开 → 更接近等概
- 低温（\(\beta\) 大）→ 约束紧 → 分布集中在低能态
- \(T \to \infty\) → \(\beta \to 0\) → \(p \to \text{Uniform}\)（所有态等概）
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")
# 左图：不同温度的能量分布
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)
# 右图：温度 vs 熵
# 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 温度是什么？\(1/T = \partial S / \partial E\)
“\(\lambda\) 就是 \(1/k_BT\)”不只是换个记号——它是热力学温度定义和最大熵乘子的严格等同。对最大熵解 \(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\]
对 \(\langle E\rangle\) 求导（注意 \(Z\) 通过 \(\lambda\) 依赖 \(\langle E\rangle\)，用 \(\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\]
而热力学里温度的定义（Clausius，\(dS = \delta Q / T\)）正是 \(\dfrac{1}{T} = \dfrac{\partial S}{\partial E}\)。两式对照（\(S = k_B H\)）：
\[\lambda = \frac{1}{k_B T}\]
用本节的指数分布验证：\(H = 1 + \ln\langle E\rangle\)，\(\frac{dH}{d\langle E\rangle} = \frac{1}{\langle E\rangle} = \frac{1}{k_BT}\)。✓
物理直觉：温度衡量”每注入一单位能量，熵涨多少”——冷系统涨得多（\(1/T\) 大），热系统涨得少。两个系统接触时，能量从热流向冷，是因为同一份能量在冷系统那边换到的熵更多，总熵增大——热传导的方向也是最大熵推出来的，不需要额外假设。
6.3 线性能级 → 指数分布
当能量态密度均匀（如量子谐振子的等间距能级，或二维气体的动能），归一化后：
\[P(E) = \beta\, e^{-\beta E} = \text{Exponential}(\beta = 1/k_BT)\]
这就是指数分布。Boltzmann 分布和指数分布不是”类似”——它们是同一个东西。
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：反应速率是 Boltzmann 尾巴
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")
化学反应的 Arrhenius 公式 \(k = A\, e^{-E_a/(k_BT)}\) 就是 Boltzmann 分布的尾部概率——能越过能垒 \(E_a\) 的分子比例。
6.4 二次动能 → 正态速度分布
三维气体分子的动能是速度的二次函数：\(E_x = \frac{1}{2}mv_x^2\)
将 Boltzmann 因子作用到动能上：
\[P(v_x) \propto e^{-mv_x^2 / (2k_BT)} = e^{-v_x^2 / (2\sigma^2)}, \quad \sigma = \sqrt{k_BT/m}\]
这正是正态分布。对比最大熵视角：\(f(x) = x^2\)，拉格朗日乘子 \(\lambda_2 = m/(2k_BT) = 1/(2\sigma^2)\)——物理里的 \(m/(2k_BT)\) 就是最大熵推导里的 \(\lambda_2\)。
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 速率分布
速率 \(v = \sqrt{v_x^2 + v_y^2 + v_z^2}\)，由三个独立正态分量合成：
\[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 配分函数：不只是归一化常数
前面把 \(Z\) 一直当成”把概率凑成 1 的那个分母”。它其实是这一章里最能干的一件工具。
\[Z(\beta) \;=\; \sum_i e^{-\beta E_i} \qquad\text{（连续情形换成积分）}\]
把 \(\ln Z\) 对 \(\beta\) 求导，平均能量自己掉出来：
\[-\frac{\partial \ln Z}{\partial \beta} = \frac{1}{Z}\sum_i E_i\,e^{-\beta E_i} = \bar E\]
再求一次导，掉出来的是涨落：
\[\frac{\partial^2 \ln Z}{\partial \beta^2} \;=\; \overline{E^2} - \bar E^2 \;=\; \mathrm{Var}(E)\]
一个函数 \(\ln Z\)，求一次导给平均值，求两次导给方差。 这不是巧合——它是矩母函数那套机器：\(\ln Z(\beta)\) 就是能量的累积量母函数。
拿最小的系统试一遍：两能级系统（能量 \(0\) 与 \(\Delta\)，比如一个只能”基态/激发态”二选一的分子）。\(Z = 1 + e^{-\beta\Delta}\)，于是
\[p_1 = \frac{e^{-\beta\Delta}}{1+e^{-\beta\Delta}}, \qquad \bar E = \Delta\,p_1\]
D <- 1 # 能隙 Δ，取 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，数值导数
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)
左图：低温时系统冻在基态（\(p_1\to0\)，没有热预算去爬那个能隙），高温时两个能级平分（\(p_1\to1/2\)，\(k_BT\) 远大于 \(\Delta\)，能隙形同不存在）。注意高温极限不是”全跑到高能级”，而是等概率——这正是 \(T\to\infty\) 时玻尔兹曼因子趋于均匀。
右图是这一节的重点。热容 \(C_V = \mathrm{d}\bar E/\mathrm{d}T\) 是你能测的宏观响应：升一度要灌多少能量。而 \(\mathrm{Var}(E)\) 是微观涨落。两条曲线严丝合缝地压在一起，因为
\[\boxed{\;\mathrm{Var}(E) \;=\; k_B T^2\, C_V\;}\]
——涨落和响应是同一个二阶导的两种读法（\(\partial^2\ln Z/\partial\beta^2\)）。系统越”容易被加热”，它的能量本来就抖得越厉害。 峰出现在 \(k_BT \approx 0.417\,\Delta\)（Schottky 峰）：温度太低爬不上去、太高已经饱和，只有在能隙尺度附近，加一点热量才最能改变占据数。
这台”\(\ln Z\) 求导出矩”的机器不是热力学专属。记住它的三个零件——指数里的 \(\beta\)、被平均的 \(E\)、归一化的 \(\ln Z\)——末章会看到它们在统计学里各有名字（自然参数、充分统计量、累积量母函数），而玻尔兹曼分布只是这套结构最早被物理学发现的那个特例。
6.7 简并与态密度：玻尔兹曼因子只是一半
到这里有个容易被忽略的缺口。玻尔兹曼因子说的是”每一个状态的概率 \(\propto e^{-\beta E}\)“。但如果你关心的不是状态而是能量，就还得数清楚：同一个能量上有多少个状态。
\[p(E) \;\propto\; \underbrace{g(E)}_{\text{有多少个格子}}\;\cdot\;\underbrace{e^{-E/k_BT}}_{\text{每个格子多重}}\]
\(g(E)\) 是简并度（离散能级）或态密度（连续能量）。两个因子的方向正好相反——格子数往上、权重往下——于是概率的峰不在能量最低处：
Tt <- 1; b <- 1/Tt; s <- 2 # g(E) ∝ E^s，取 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)
把 \(g(E)\propto E^s\) 代进去求极值，峰的位置有闭式：
\[\frac{\mathrm{d}}{\mathrm{d}E}\Bigl[s\ln E - \beta E\Bigr] = 0 \;\Longrightarrow\; E^\ast = \frac{s}{\beta} = s\,k_BT\]
图里 \(s=2\)、\(k_BT=1\)，峰精确落在 \(E^\ast=2\)。所以”最可能的能量”随温度线性上移，而且恒不为零——尽管每一个具体状态都是越低越可能。这就是本章前面 Maxwell 速率分布那个 \(v^2\) 因子的来处：三维速度空间里半径为 \(v\) 的球壳面积 \(\propto v^2\)，就是 \(g(v)\)。速度分量是纯高斯（\(g\) 均匀），速率却有个先升后降的峰。
这里请记住这个 \(g\) 的位置：它不来自约束、也不来自玻尔兹曼因子，而是”你拿什么当计数单位”。末章会给它一个统计学的名字（基测度），并说明它其实一直潜伏在前面每一个推导里。
6.8 等分定理：每个二次自由度分到 \(\frac12 k_BT\)
第5章说过：能量对坐标二次，分布就是高斯。反过来问一句——那这个自由度平均拿到多少能量?
一维弹簧势 \(E = \frac12\kappa x^2\)，位置服从 \(\sigma^2 = k_BT/\kappa\) 的高斯，于是
\[\bar E = \tfrac12\kappa\,\overline{x^2} = \tfrac12\kappa\sigma^2 = \tfrac12 k_BT\]
\(\kappa\) 消掉了。 弹簧软还是硬、分子轻还是重，都不影响这个自由度分到的能量——只由温度决定。
kap <- 2.7 # 任取一个刚度，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推广就是等分定理：能量里每一个二次项（每个方向的动能、每个简谐坐标）平均分到 \(\frac12 k_BT\)；\(f\) 个二次自由度给出 \(\bar E = \frac f2 k_BT\)、\(C_V = \frac f2 k_B\)。单原子气体三个平动方向 → \(\frac32 k_BT\)，这就是理想气体那个熟悉的结果。
顺带解释了一件事：为什么温度可以用”每个自由度的能量”来定义。\(k_BT\) 不是随便一个换算因子，它就是热浴发给每个二次自由度的那份”零花钱”。
6.9 自由能：最大熵的另一面
最后补一个抽象台阶——它是本章通向末章的桥。
把 \(\bar E = -\partial_\beta\ln Z\) 和熵合起来算（\(S = k_B H\)，\(H = -\sum p\ln p\)，代入 \(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\) 就是自由能。于是同一件事有了两种说法：
- 最大熵视角：固定平均能量 \(\bar E\)，最大化熵 \(S\)；
- 自由能视角：固定温度 \(T\)，最小化 \(F = \bar E - TS\)。
两者是一对 Legendre 变换（\(\bar E\) 与 \(\beta\) 互为共轭变量，正如价格与数量）。这解释了物理里那句常见的话”系统趋于最小自由能”——它不是与最大熵原理并列的另一条定律，而是同一条定律在”温度给定”而非”能量给定”时的写法。
\(F\) 里那个减号也就有了含义：\(-TS\) 是熵得到的补贴。温度越高，多样性越值钱；\(T\to0\) 时熵项失效，系统只认能量最低。这正是第1章”温度衡量每注入一单位能量、熵涨多少”的对偶面。
末章会看到 \(\ln Z\) 在统计学里的名字，以及这个 Legendre 对偶如何变成指数族的凸共轭——同一套结构，两个学科各自命名了一遍。
7 统一图景：指数族与 \(p \propto e^{-\lambda f(x)}\)
7.1 一个方程，四种分布
回到最大熵的通解 \(p \propto e^{-\lambda f(x)}\)：
\[p(x) \propto e^{-\lambda_1 f_1(x) - \lambda_2 f_2(x) - \cdots}\]
这在统计学中称为指数族（exponential family）。自然界反复出现这些分布，不是因为它们”好看”，而是因为它们是在有限信息下最不偏不倚的选择。
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 先说清那把尺子：测度是什么
下一节要用到”参照测度”。这个词听着像泛函分析的门槛，其实东西朴素得多，而且不讲清它，下一节的主角就会变成一句咒语。先花一节把它说透。
测度就是”你怎么数有多少”。 给你一段区间 \([a,b]\)，问它”有多大”，最自然的回答是长度 \(b-a\)。把这种”给集合指定大小”的规则严格化的人是 Henri Lebesgue（1902），所以”按普通长度来数”这把尺子就叫 Lebesgue 测度，在积分里就是那个不起眼的 \(\mathrm{d}x\)。
关键在于：尺子不止一把。 你可以按长度数，也可以按”数量级”数（每十倍算一格），也可以在整数上一个点算一格。换一把尺子，“均匀”“分散”“无知”这些词的含义就跟着变。
7.2.1 为什么熵必须先说清尺子
离散熵不用操心这件事：它数的是状态，无量纲，换什么单位都一样。连续情形不然——微分熵不是坐标无关的。看一个例子：
Hnum <- function(dens, lo, hi) # 数值积分算微分熵
-integrate(function(x){p <- dens(x); ifelse(p > 0, p*log(p), 0)}, lo, hi)$value
# 同一根"0 到 1 米之间均匀随机"的长度，只是换记录单位
c(以米计 = Hnum(function(x) dunif(x, 0, 1), 0, 1),
以厘米计 = Hnum(function(x) dunif(x, 0, 100), 0, 100),
以毫米计 = Hnum(function(x) dunif(x, 0, 1000), 0, 1000),
ln100 = log(100), ln1000 = log(1000)) |> round(4)## 以米计 以厘米计 以毫米计 ln100 ln1000
## 0.0000 4.6052 6.9078 4.6052 6.9078同一件随机事实，换个单位，熵的数值就变了——以米计是 \(0\)，以厘米计变成 \(\ln 100 = 4.605\)。原因不神秘：概率密度 \(p\) 是带单位的（每米多少概率），换单位时 \(p\) 整体缩放，\(\ln p\) 就整体平移。
所以”微分熵的绝对数值”本身没有意义，它只在给定一把尺子之后才有意义。这也顺带解释了前面见过的那件怪事：指数分布 \(\mu < 1/e\) 时 \(H<0\)——负号不代表”信息量为负”，只代表这个分布比”单位长度那么宽”更集中。微分熵测的是”相对于尺子有多分散”。
7.2.2 相对熵才是不依赖尺子的那个量
既然如此，正确的做法是把尺子写进公式。相对熵（KL 散度）就是这么干的：
\[D(p\,\|\,m) \;=\; \int p(x)\,\ln\frac{p(x)}{m(x)}\,\mathrm{d}x\]
\(m\) 就是那把尺子（参照测度）。换单位时 \(p\) 和 \(m\) 同步缩放，比值 \(p/m\) 不变，于是 \(D\) 也不变：
# p = Beta(2,5)；参照 m = 该区间上的均匀分布。再换成厘米看一遍
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(相对熵_米 = D_m, 相对熵_厘米 = D_cm,
微分熵_米 = Hnum(function(x) dbeta(x,2,5), 0, 1),
微分熵_厘米 = Hnum(function(y) dbeta(y/100,2,5)/100, 0, 100)) |> round(5)## 相对熵_米 相对熵_厘米 微分熵_米 微分熵_厘米
## 0.48453 0.48453 -0.48453 4.12064相对熵两次完全相同，微分熵却差了 \(\ln 100\)。所以严格地说，该最大化的从来是相对熵；本文前面一路写的”最大化 \(H\)“，是”参照测度取成常数（Lebesgue）“时的简写。 这句话就是下一节整个推广的全部内容。
7.2.3 术语对表（同一个东西的几个名字）
这里有个纯粹的命名混乱，值得单独澄清一下：
| 说法 | 指什么 |
|---|---|
| 参照测度 / 背景测度 / 基测度 / base measure | 都是同一个东西：公式里那个 \(m\)，你选的那把尺子 |
| Lebesgue 测度 | 尺子的一个具体选择：按普通长度数（\(m=\) 常数，积分里的 \(\mathrm{d}x\)） |
| 计数测度 | 另一个具体选择：离散点上，一个点算一格（\(m(y)=1\)） |
| 态密度 \(g(E)\) / 简并度 | 物理学对同一件事的叫法（第 6 章那个 \(g\)） |
也就是说：“参照测度”是角色，“Lebesgue 测度”是演员之一。
7.2.4 三把常用的尺子
| 尺子 | 数法 | “什么都不知道”时给出 |
|---|---|---|
| 等长尺（Lebesgue，\(\mathrm{d}x\)） | 每段等长同等重要，普通直尺 | 有界区间上的均匀分布 |
| 计数尺（\(m(y)=1\)） | 整数上一个点算一格 | 整数上的均匀；配固定均值给几何分布 |
| 对数尺（\(\mathrm{d}x/x\)） | 每个数量级等宽，像地震震级、pH、分贝 | \(p\propto 1/x\)，一条幂律 |
第三把最值得试一下手感。\(\mathrm{d}x/x\) 是”对数刻度上的等长”（因为 \(\mathrm{d}\ln x = \mathrm{d}x/x\)）。在这把尺子上均匀，意思是每个数量级的概率相同：
Z <- log(1000) # 支撑取 [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三个数量级各占 \(1/3\)——落在 1 到 10 之间和落在 100 到 1000 之间一样可能。这正是”对数尺度上什么都不知道”的样子，而它在普通直尺上看就是幂律 \(1/x\)。为什么这把尺子配幂律天然合适？因为 \(\mathrm{d}x/x\) 是唯一在 \(x\to cx\)（整体放大）下不变的尺子——量纲上”倍数”才是它的单位，而幂律正是没有特征尺度的那个分布。
准备工作到此为止。下一节把 \(m\) 正式请进公式，看它如何一口气点亮整个分布动物园。
7.3 完整的指数族：一个模板，点亮整个动物园
前面每个分布都写成 \(p\propto e^{-\lambda f(x)}\)。但这还不是最全的形式——它悄悄假设了一样东西。补上它，就得到统计学里统一万物的那个对象：指数族。
被隐藏的参照测度。 写 \(H=-\int p\ln p\,dx\) 时，其实默认了一个参照：\(m(x)=\) 常数（Lebesgue 测度，“每段等长的 \(dx\) 同等重要”）。更一般地，该最大化的是相对熵 \(-\int p\ln\frac{p}{m}\,dx\)；对 \(p\) 逐点求变分（\(\frac{\delta}{\delta p}\bigl[-p\ln\frac{p}{m}\bigr] = -\ln\frac{p}{m}-1\)，令它加上乘子项为零）解出
\[\boxed{\,p(x)\;\propto\;m(x)\,e^{-\sum_i \lambda_i f_i(x)}\,}\]
这就是指数族的一般形式。两个零件生成整个分布动物园：参照测度 \(m\)（你在什么”背景几何”上计数）＋ 约束函数 \(\{f_i\}\)（你知道哪些量）。前面各章都是它取 \(m=\) 常数、\(f\) 取 \(x\)／\(x^2\)／\(\ln x\)／\(E\) 的特例。
换一个背景测度，就换出一整族分布。
离散计数 \(y=0,1,2,\dots\)：最朴素的情形——非负整数 + 固定均值，普通最大熵给的不是 Poisson，而是 Geometric \(p_y=(1-\theta)\theta^y\)（指数分布在整数格点上的孪生）。要得到 Poisson，必须把参照测度换成 \(m(y)=1/y!\)；Binomial 要 \(m(y)=\binom{n}{y}\)。这些组合数是”计数的背景结构”，约束变不出来，只能由 \(m\) 提供：
| 分布 | 基测度 \(m(y)\) | 约束 | 结果 |
|---|---|---|---|
| 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}\) |
尺度不变测度 \(dx/x\)（对 \(\ln x\) 均匀，因为 \(d\ln x=dx/x\)）：什么约束都不加，解就是 \(p\propto 1/x\)——幂律。所以幂律是”对数尺度上什么都不知道”的分布，正如均匀是”线性尺度上什么都不知道”。base measure 不是补丁，而是最大熵一直在用的隐藏坐标系：换掉它，“最无知”长什么样就跟着变。
格局。 \(m(x)\,e^{-\sum\lambda_i f_i}\) 不是本文的小把戏，它是现代统计的核心对象——指数族：\(\{f_i\}\) 是充分统计量，\(\lambda_i\) 是自然参数；广义线性模型（GLM）、共轭先验、最大似然的一整套漂亮性质都建在这个形式上。回头看，本文一路点亮的指数、幂律、正态、玻尔兹曼、几何、泊松……不过是同一个母体——指数族——的不同角落。
7.4 为什么自然界总是这些分布？
不是因为自然”选择”了这些函数形式，而是因为：
在约束条件下，微观实现方式最多的宏观状态就是最大熵状态。
自然界不追求最低能量，也不追求最高能量，而是追求最容易实现的状态——即微观排列数最多的分布。指数型衰减 \(e^{-\lambda f(x)}\) 之所以反复出现，是因为它是约束优化问题的通解。
8 全文总结
一条主线贯穿全篇：优化一个量 + 施加若干约束 → 唯一确定的结构。
- 优化什么：熵 \(H = -\int p\ln p\,dx\)。它不是任意选的——计数论证证明了 \(\ln W \approx N\cdot H\)，熵最大的分布就是微观排列数压倒性最多、因而几乎必然出现的那个（\(W \approx e^{NH}\)）。
- 怎么优化：拉格朗日乘子 / 变分法。约束吃进拉格朗日函数，对分布逐点求导，通解一步落地：\(p(x) \propto e^{-\lambda f(x)}\)。
- 约束是什么：由你的已知信息（测量或守恒律）决定。换一个约束函数 \(f(x)\)，就换一个分布：
| 约束 \(f(x)\) | 得到的分布 | 乘子的物理身份 | 实例 |
|---|---|---|---|
| 无（仅有界） | 均匀分布 | — | — |
| \(x\)（均值，非负） | 指数分布 | — | 等温大气、放射性衰变 |
| \(\ln x\)（几何均值，\(x\ge x_{\min}\)） | 幂律分布 / Pareto | — | 城市规模、财富尾部、网络度数 |
| \(x,\ x^2\)（均值+方差） | 正态分布 | — | 化学键振动、布朗扩散 |
| \(E(x)\)（能量） | 玻尔兹曼分布 | \(\lambda = 1/k_BT\)（温度） | Maxwell 速率、Arrhenius |
指数、幂律、正态、玻尔兹曼——这四个看似无关的分布，因此是同一个方程在不同约束下的四张脸（均匀分布则是”无约束”的基线）。而这套”优化 + 约束”的模板远不止于此：拉格朗日乘子在力学里是约束力、在热力学里是温度/化学势/压强、在优化里是对偶变量——它是连接数学、物理、化学的一根共同轴。
一句话带走：指数族不是自然界的审美偏好，而是”给定参照测度、只知道有限几个平均量”这一推断处境的数学结论。
这句话的分量要放准：它的必然性属于推断，不属于自然。三个前提（参照测度固定、约束有限且已知、最大化熵）任一松动，结论就不成立——自然界里确有不属于指数族的分布：Student-\(t\)、Cauchy、Lévy 稳定律、混合分布、以及支撑随参数变动的族（未知 \(x_{\min}\) 的 Pareto、未知形状的 Weibull）。反过来也一样要紧：是指数族的分布，也常常不是靠”诚实推断”产生的——正态可以由中心极限定理（加法聚合）长出来，幂律可以由优先增长、Kesten 过程或临界分枝长出来，这些机制与”在有限信息下最诚实”没有关系。
所以同一个分布往往有两套互不隶属的来路：最大熵回答”在这点信息下我该假设什么”，机制回答”它实际会长成什么”。 两条路频繁重合（正态、指数、幂律都各有两套解释），但它们回答的不是同一个问题——重合本身才是这件事最耐琢磨的地方。
作者：胡杨博士 ｜ 首次发布：2026-07-12 ｜ 最后更新：2026-08-26 16:24 EDT
参考文献
留言与讨论
读到不清楚、觉得有错、或想到别的例子，都欢迎在下面留言（用 GitHub 账号登录即可）。
留言支持 Markdown、代码块和 $LaTeX$ 公式。