第 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\) 时对应的税额、补贴额或福利指标。只要分布稍微复杂，这些积分就未必有解析解。计算统计要回答的是：当积分不能手工算出时，怎样用可控误差的数值方法近似它？

本章有两个重点。第一，理解一维数值积分如何把“连续面积”近似为有限个函数值的加权和；第二，理解 Laplace 近似如何利用最大点附近的局部二次形状近似一个整体积分。前者适合低维、函数值容易计算的问题，后者常用于似然、后验分布和边际似然的近似计算。

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\) 越大，近似通常越精细，但计算量也越大。

为了看清数值积分的计算结构，先手写一个梯形公式函数。它只适用于有限区间，但足以说明思想。

trapezoid_integral <- function(f, a, b, m = 1000) {
  x <- seq(a, b, length.out = m + 1)
  y <- f(x)
  delta <- (b - a) / m
  delta * (sum(y) - 0.5 * y[1] - 0.5 * y[m + 1])
}

下面用一个模拟的收入分布例子说明。假设家庭收入 \(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

# 手写梯形公式
p_trapezoid <- trapezoid_integral(d_income, a = 0, b = poverty_line, m = 2000)

c(exact = p_exact,
  integrate = p_integrate,
  trapezoid = p_trapezoid)
#>     exact integrate trapezoid 
#>     0.124     0.124     0.124

R 中的 integrate() 会自动调整积分网格，比手写固定网格更稳健。对本科阶段的计算统计而言，需要掌握两点：第一，积分在统计推断中经常出现；第二，当有可靠的标准函数时，应优先使用它，而不是盲目从零实现复杂数值积分。

网格越密，梯形公式通常越准确。下面比较不同网格数量下的误差。

m_grid <- c(20, 50, 100, 500, 2000)
trap_error <- sapply(m_grid, function(m) {
  abs(trapezoid_integral(d_income, 0, poverty_line, m) - p_exact)
})
data.frame(m = m_grid, absolute_error = trap_error)
#>      m absolute_error
#> 1   20      6.577e-05
#> 2   50      1.052e-05
#> 3  100      2.631e-06
#> 4  500      1.052e-07
#> 5 2000      6.577e-09

这个表说明了一个重要事实：数值积分不是“调用一次函数就结束”，而是要考虑区间、网格、误差和计算量。在统计应用中，如果被积函数有尖峰、厚尾或不连续点，固定网格可能会表现很差，此时更应使用自适应积分、变量变换或随机模拟方法。

如果要计算一个非线性政策函数的期望，也可以直接把函数写进被积函数中。下面假设收入低于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 近似常出现为

\[ I_n=\int a(x)\exp\{n h(x)\}\,dx, \]

其中 \(n\) 往往对应样本量，\(h(x)\) 是平均对数似然或平均对数后验的一部分，\(a(x)\) 是变化相对平缓的函数。如果 \(h(x)\) 在 \(x_0\) 处达到内部最大值，则

\[ I_n\approx a(x_0)\exp\{n h(x_0)\} \sqrt{\frac{2\pi}{n[-h''(x_0)]}}. \]

这个形式特别有助于理解为什么样本量越大，Laplace 近似往往越好：当 \(n\) 增大时，\(\exp\{n h(x)\}\) 会更集中在最大点附近，远离最大点的区域对积分贡献迅速变小。

一维 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)\)；
代入近似公式。

实际编程时，建议在对数尺度上实现 Laplace 近似。下面写一个一维辅助函数：给定log_g(x)，先用 optimize() 找最大点，再用有限差分近似二阶导数。

laplace_1d <- function(log_g, lower, upper, eps = 1e-4) {
  opt <- optimize(log_g, interval = c(lower, upper), maximum = TRUE)
  x0 <- opt$maximum
  h0 <- opt$objective
  h2 <- (log_g(x0 + eps) - 2 * log_g(x0) + log_g(x0 - eps)) / eps^2

  if (!is.finite(h2) || h2 >= 0) {
    stop("最大点附近的二阶导数必须为负。")
  }

  list(
    mode = x0,
    second_derivative = h2,
    log_integral = h0 + 0.5 * log(2 * pi) - 0.5 * log(-h2)
  )
}

这个函数只适合最大点不在边界附近、且 lower 到 upper 覆盖主要概率质量的情形。它的作用是帮助理解算法步骤；正式项目中仍要结合图形、多个区间和数值检查。

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,
  log_error = log_int_laplace - log_int_exact)
#>                 mode           exact_mean laplace_log_integral 
#>               0.3113               0.3122            -255.8643 
#>   exact_log_integral            log_error 
#>            -255.8660               0.0017

在这个样本量下，Laplace 近似给出的归一化常数已经比较接近解析结果。需要注意的是，Laplace 近似使用的是峰值附近的信息，因此它对单峰、近似对称、样本量较大的后验分布更可靠。如果后验分布偏斜很强、峰值靠近边界，或者存在多个峰，近似误差可能较大。

对同一个例子，也可以把后验分布近似为正态分布。由于

\[ h''(\theta_0)<0, \]

后验在 \(\theta_0\) 附近的近似方差为

\[ \operatorname{Var}(\theta\mid x)\approx -\frac{1}{h''(\theta_0)}. \]

下面把 Laplace 近似区间与精确 Beta 后验区间比较。

theta_sd_laplace <- sqrt(-1 / h2(theta_mode))

c(laplace_q025 = theta_mode + qnorm(0.025) * theta_sd_laplace,
  laplace_q975 = theta_mode + qnorm(0.975) * theta_sd_laplace,
  exact_q025 = qbeta(0.025, a_post, b_post),
  exact_q975 = qbeta(0.975, a_post, b_post))
#> laplace_q025 laplace_q975   exact_q025   exact_q975 
#>       0.2663       0.3562       0.2683       0.3578

这里两组区间相近，是因为样本量不小且后验峰值远离 0 和 1。如果样本量很小或样本比例接近 0、1，直接在 \(\theta\) 尺度上做正态近似就可能给出小于 0 或大于 1 的区间端点。此时更稳妥的做法是对logit 变换后的参数做近似，或者使用网格、Monte Carlo、MCMC 等方法。

7.5 案例：企业订单率的边际似然近似

下面看一个没有简单解析解的例子。假设我们观察到某类小微企业在若干月份的线上订单数\(y_1,\ldots,y_n\)，并用 Poisson 模型描述：

\[ y_i\mid \eta \sim \operatorname{Poisson}\{\exp(\eta)\}, \]

其中 \(\eta=\log \lambda\)，\(\lambda\) 是平均订单数。为了表达对订单率的先验认识，设

\[ \eta\sim N(m_0,s_0^2). \]

边际似然为

\[ p(y)= \int \left[ \prod_{i=1}^n \frac{\exp\{-\exp(\eta)\}\exp(\eta y_i)}{y_i!} \right] \phi(\eta;m_0,s_0^2) \,d\eta. \]

其中 \(\phi(\eta;m_0,s_0^2)\) 是均值为 \(m_0\)、方差为 \(s_0^2\) 的正态密度。这个积分通常没有闭式解，但被积函数往往在后验众数附近高度集中，适合使用 Laplace 近似。

set.seed(2027)
y_orders <- rpois(40, lambda = 3.5)
m0 <- log(3)
s0 <- 0.7

log_joint_eta <- function(eta) {
  sum(dpois(y_orders, lambda = exp(eta), log = TRUE)) +
    dnorm(eta, mean = m0, sd = s0, log = TRUE)
}

lap <- laplace_1d(log_joint_eta, lower = -1, upper = 3)

integrand_eta <- function(eta) {
  sapply(eta, function(e) exp(log_joint_eta(e)))
}
log_marginal_numeric <- log(integrate(integrand_eta,
                                      lower = -Inf,
                                      upper = Inf)$value)

c(mode_eta = lap$mode,
  mode_lambda = exp(lap$mode),
  log_marginal_laplace = lap$log_integral,
  log_marginal_numeric = log_marginal_numeric,
  log_error = lap$log_integral - log_marginal_numeric)
#>             mode_eta          mode_lambda log_marginal_laplace 
#>              1.30548              3.68946            -81.44543 
#> log_marginal_numeric            log_error 
#>            -81.41671             -0.02872

在这个例子中，数值积分仍然可行，因此我们可以用它检查 Laplace 近似。若参数维度更高，直接数值积分会迅速变得困难，Laplace 近似就能提供一个相对便宜的近似方案。

这个案例也说明了第 6 章和本章的联系：Laplace 近似的第一步是寻找对数被积函数的最大点，这本质上是一个优化问题；第二步是计算最大点附近的曲率，这又依赖 Hessian 或数值二阶导数。

7.6 多维情形与维数困难

如果参数是 \(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。

7.7 实践中的注意事项

数值积分和 Laplace 近似在教材中看起来整洁，但实际使用时要注意以下问题。

优先在对数尺度上计算。似然函数和后验密度可能非常小，直接计算容易下溢。先计算\(\log g(x)\)，最后再还原，通常更稳定。
检查积分区间。如果主要概率质量不在所选区间内，任何网格积分或优化都会给出误导结果。
注意边界最大点。Laplace 近似的标准公式假设最大点在内部；若最大点靠近边界，正态近似可能很差。
警惕多峰分布。一个局部最大点不能代表整个积分。必要时应画图或尝试多个初始区间。
区分算法误差和统计误差。数值积分误差来自计算近似；抽样误差来自数据随机性，二者含义不同。
用可检查例子校准直觉。像 Beta–Binomial 这类有解析解的例子很适合用来检查近似方法是否合理。

7.8 本章小结

本章介绍了统计计算中的积分问题。一维数值积分可以用网格、梯形公式或 R 中的 integrate() 近似；Laplace 近似则把积分主要贡献集中在峰值附近，用二次展开和正态形状近似积分。它在极大似然近似、贝叶斯后验近似和边际似然计算中都很有用。Beta–Binomial 例子帮助我们用解析结果检查近似误差，Poisson 订单率案例展示了 Laplace 近似如何处理没有简单闭式解的边际似然。与此同时，本章也提醒我们：积分维度升高后，普通网格法很快不可行，这正是后续 Monte Carlo 和 MCMC 方法的重要动机。

7.9 练习

令 \(Y\) 服从对数正态分布，改变 \(\mu\) 和 \(\sigma\)，用 integrate() 计算 \(P(Y<c)\)，并解释参数变化对低收入比例的影响。
对函数 \(h(x)=-x^4/4+x^2/2\)，画出 \(\exp\{h(x)\}\) 的图形，并讨论 Laplace 近似是否合适。
在 Beta–Binomial 例子中，把样本量从 40 增加到 4000，比较 Laplace 近似误差的变化。
在 Poisson 订单率案例中，把先验标准差 \(s_0\) 改为 0.2 和 2，比较 Laplace 近似误差。
思考为什么二维网格积分比一维积分困难得多。设每个维度取 1000 个点，二维和五维分别需要多少个函数值？

第 7 章 数值积分与近似方法