第 7 章 数值积分与近似方法
上一章讨论了如何通过优化寻找目标函数的最大值或最小值。本章转向另一类同样常见的计算问题：积分。在统计学中，积分并不只是数学分析里的练习题，而是期望、概率、边际似然、后验均值和预测分布的共同语言。
例如，在收入分布研究中，若 \(Y\) 表示家庭年收入，\(f(y)\) 是其概率密度函数，我们可能关心低收入比例
\[ P(Y<c)=\int_0^c f(y)\,dy, \]
也可能关心某个税收或补贴规则 \(g(Y)\) 的平均财政影响
\[ E\{g(Y)\}=\int_0^\infty g(y) f(y)\,dy. \]
这里 \(c\) 是低收入线，\(g(y)\) 是收入为 \(y\) 时对应的税额、补贴额或福利指标。只要分布稍微复杂，这些积分就未必有解析解。计算统计要回答的是：当积分不能手工算出时，怎样用可控误差的数值方法近似它？
7.1 积分在统计推断中的作用
设 \(X\) 是随机变量，概率密度为 \(f(x)\)。函数 \(g(X)\) 的期望为
\[ \mathbb{E}\{g(X)\}=\int g(x)f(x)\,dx. \]
其中，\(x\) 表示随机变量可能取到的值，\(f(x)\) 表示这些取值的相对可能性，\(g(x)\) 表示我们真正关心的统计量或经济指标。许多统计问题都可以写成这种形式：
- 分布函数 \(F(c)=P(X\leq c)=\int_{-\infty}^c f(x)\,dx\)；
- 后验均值 \(E(\theta\mid y)=\int \theta p(\theta\mid y)\,d\theta\)；
- 预测分布 \(p(\tilde y\mid y)=\int p(\tilde y\mid\theta)p(\theta\mid y)\,d\theta\)；
- 边际似然 \(p(y)=\int p(y\mid\theta)p(\theta)\,d\theta\)。
这些式子中的积分有时可以解析求解。例如正态分布、Beta–Binomial 模型、Normal–Normal 模型等有漂亮的闭式结果。但在多数实际模型中，积分需要借助数值方法、Monte Carlo 方法或近似方法完成。本章先介绍确定性的数值积分和 Laplace 近似，后面第 10 章会介绍 Monte Carlo 积分。
7.2 一维数值积分
一维积分最直观的想法是把积分区间切成许多小段，在每一段上用简单图形近似曲线下面积。以区间\([a,b]\) 上的积分为例：
\[ I=\int_a^b h(x)\,dx. \]
若把区间分成 \(m\) 个等长小区间，步长为 \(\Delta=(b-a)/m\)，网格点为\(x_j=a+j\Delta\)，梯形公式为
\[ I \approx \Delta\left\{ \frac{h(x_0)}{2}+\sum_{j=1}^{m-1}h(x_j)+\frac{h(x_m)}{2} \right\}. \]
这里 \(h(x)\) 是被积函数，\(m\) 是区间划分数量。\(m\) 越大，近似通常越精细，但计算量也越大。
下面用一个模拟的收入分布例子说明。假设家庭收入 \(Y\) 服从对数正态分布，即 \(\log Y\sim N(\mu,\sigma^2)\)。这不是某个真实数据来源，而是一个用于教学的模拟设定。
mu <- log(8)
sigma <- 0.6
poverty_line <- 4
d_income <- function(y) dlnorm(y, meanlog = mu, sdlog = sigma)
# 低收入比例：解析分布函数
p_exact <- plnorm(poverty_line, meanlog = mu, sdlog = sigma)
# 数值积分
p_integrate <- integrate(d_income, lower = 0, upper = poverty_line)$value
c(exact = p_exact, numerical = p_integrate)
#> exact numerical
#> 0.124 0.124R 中的 integrate() 会自动调整积分网格，比手写固定网格更稳健。对本科阶段的计算统计而言，需要掌握两点：第一，积分在统计推断中经常出现；第二，当有可靠的标准函数时，应优先使用它，而不是盲目从零实现复杂数值积分。
如果要计算一个非线性政策函数的期望，也可以直接把函数写进被积函数中。下面假设收入低于4 万元的家庭获得补贴，补贴金额随收入线性递减：
\[ g(y)= \begin{cases} 2(1-y/4), & 0<y<4,\\ 0, & y\geq 4. \end{cases} \]
其中 \(g(y)\) 的单位可以理解为万元。
subsidy <- function(y) ifelse(y < poverty_line, 2 * (1 - y / poverty_line), 0)
cost <- integrate(function(y) subsidy(y) * d_income(y),
lower = 0, upper = Inf)$value
cost
#> [1] 0.05831这个数值就是每户平均补贴支出的理论近似。实际调查中，\(f(y)\) 未必已知，后面我们会用 Monte Carlo和 Bootstrap 处理类似问题。
7.3 Laplace 近似的基本思想
许多统计积分具有如下形式：
\[ I=\int g(x)\,dx, \]
其中 \(g(x)>0\)，并且大部分面积集中在某个峰值附近。Laplace 近似的思想是：如果 \(g(x)\) 在最大点\(x_0\) 附近最重要，那么可以只用 \(x_0\) 附近的二次近似来估计整个积分。
令
\[ h(x)=\log g(x). \]
如果 \(g(x)\) 在 \(x_0\) 处取得最大值，那么 \(h(x)\) 也在 \(x_0\) 处取得最大值，并且\(h'(x_0)=0\)。在 \(x_0\) 附近做二阶泰勒展开：
\[ h(x)\approx h(x_0)+\frac{1}{2}h''(x_0)(x-x_0)^2. \]
由于 \(x_0\) 是局部最大点，通常有 \(h''(x_0)<0\)。于是
\[ \begin{aligned} \int g(x)\,dx &=\int \exp\{h(x)\}\,dx\\ &\approx \exp\{h(x_0)\} \int \exp\left\{\frac{1}{2}h''(x_0)(x-x_0)^2\right\}\,dx\\ &= \exp\{h(x_0)\}\sqrt{\frac{2\pi}{-h''(x_0)}}. \end{aligned} \]
这说明 Laplace 近似本质上是在峰值附近用一个正态形状近似 \(g(x)\)。其中 \(x_0\) 决定位置，\(-1/h''(x_0)\) 决定局部宽度。
一维 Laplace 近似
设 \(g(x)>0\)，\(h(x)=\log g(x)\)。若 \(h(x)\) 在内部点 \(x_0\) 处达到最大值，并且\(h''(x_0)<0\)，则
\[ \int g(x)\,dx \approx \exp\{h(x_0)\}\sqrt{\frac{2\pi}{-h''(x_0)}}. \]
计算步骤为：
- 写出 \(h(x)=\log g(x)\)；
- 用优化算法求最大点 \(x_0\)；
- 计算二阶导数 \(h''(x_0)\)；
- 代入近似公式。
7.4 例：调查比例的后验积分
假设在一个模拟的消费意愿调查中，\(n=400\) 名受访者中有 \(x=126\) 人表示未来一年可能购买某类耐用品。令 \(\theta\) 表示总体购买意愿比例，采用 Beta 先验
\[ \theta\sim \operatorname{Beta}(a,b). \]
这里 \(a\) 和 \(b\) 可以理解为先验中的“成功”和“失败”信息量。若观测到 \(x\) 个成功、\(n-x\) 个失败，则未归一化的后验密度为
\[ g(\theta)=\theta^{a+x-1}(1-\theta)^{b+n-x-1}, \quad 0<\theta<1. \]
归一化常数是积分
\[ \int_0^1 \theta^{a+x-1}(1-\theta)^{b+n-x-1}\,d\theta. \]
这个例子有解析解，因此适合用来检查 Laplace 近似。
n <- 400
x <- 126
a <- 2
b <- 8
a_post <- a + x
b_post <- b + n - x
h <- function(theta) {
(a_post - 1) * log(theta) + (b_post - 1) * log(1 - theta)
}
h2 <- function(theta) {
-(a_post - 1) / theta^2 - (b_post - 1) / (1 - theta)^2
}
theta_mode <- (a_post - 1) / (a_post + b_post - 2)
log_int_laplace <- h(theta_mode) + 0.5 * log(2 * pi) -
0.5 * log(-h2(theta_mode))
log_int_exact <- lbeta(a_post, b_post)
c(mode = theta_mode,
exact_mean = a_post / (a_post + b_post),
laplace_log_integral = log_int_laplace,
exact_log_integral = log_int_exact)
#> mode exact_mean laplace_log_integral
#> 0.3113 0.3122 -255.8643
#> exact_log_integral
#> -255.8660在这个样本量下，Laplace 近似给出的归一化常数已经比较接近解析结果。需要注意的是，Laplace 近似使用的是峰值附近的信息，因此它对单峰、近似对称、样本量较大的后验分布更可靠。如果后验分布偏斜很强、峰值靠近边界，或者存在多个峰，近似误差可能较大。
7.5 多维情形与维数困难
如果参数是 \(d\) 维向量 \(\boldsymbol{\theta}\)，Laplace 近似可以写成
\[ \int \exp\{h(\boldsymbol{\theta})\}\,d\boldsymbol{\theta} \approx \exp\{h(\hat{\boldsymbol{\theta}})\} (2\pi)^{d/2} |\mathbf{H}|^{-1/2}, \]
其中 \(\hat{\boldsymbol{\theta}}\) 是 \(h(\boldsymbol{\theta})\) 的最大点，\(\mathbf{H}\) 是 \(-h(\boldsymbol{\theta})\) 在最大点处的 Hessian 矩阵，也就是
\[ \mathbf{H} = -\left. \frac{\partial^2 h(\boldsymbol{\theta})} {\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^\top} \right|_{\boldsymbol{\theta}=\hat{\boldsymbol{\theta}}}. \]
矩阵 \(\mathbf{H}\) 描述最大点附近后验曲面的弯曲程度。若某个方向曲率大，说明该方向的不确定性小；若曲率小，说明该方向的不确定性大。这与第 6 章中 Hessian 在优化中的作用是一致的。
但多维积分也带来“维数困难”。如果每个维度使用 \(m\) 个网格点，\(d\) 维网格就需要 \(m^d\) 个点。例如，\(m=100\)、\(d=5\) 时已经需要 \(10^{10}\) 个点，普通计算机难以承受。因此，高维积分通常不依赖朴素网格，而是使用 Laplace 近似、Monte Carlo、重要性抽样或 MCMC。