第 1 章 计算统计的动机与基本问题

在概率论、数理统计或计量经济学课程中，我们已经见过许多重要的统计公式，例如密度函数、似然函数、最小二乘估计、极大似然估计、后验分布、期望和积分。这些公式刻画了统计推断的目标，但要把它们真正用于数据分析，还必须回答一系列计算问题：公式中的量如何在计算机上求得？当矩阵规模很大、积分没有解析表达式或者估计方程不能直接求解时，统计推断应当如何继续？为什么有时公式在数学上完全正确，程序却无法运行，或者虽然正常运行，却给出不可信的结果？

这些问题构成了计算统计的研究内容。本书所称“计算统计”，是研究如何利用数值算法、随机模拟及其程序实现完成统计建模与推断，并对计算结果进行诊断和复现的学科领域。计算统计并不是用程序替代统计理论，而是研究如何把统计目标转化为可靠、稳定、高效且可验证的计算过程。简言之，统计理论说明“要计算什么”，计算统计研究“怎样可靠地把它计算出来”。在需要强调某一具体统计任务的计算实施时，本书使用“统计计算”一词。关于计算统计的基本问题和方法体系，可参见 Gentle (2009)、Monahan (2011) 和 Roger D. Peng (2018)。

本章从若干基本例子出发，说明统计推断为什么离不开计算。后续章节中的数值线性代数、求根、优化、数值积分、随机模拟、重采样和 MCMC 方法，都可以看作对本章问题的系统回答。这些方法既服务于经典统计推断，也构成贝叶斯后验计算、模型比较和后验预测的基础。

1.1 本章学习目标

完成本章后，读者应当能够：

  1. 区分统计问题、推断目标、统计模型与计算算法各自承担的作用；
  2. 说明有限精度、病态矩阵和缺少解析解为何会形成计算困难；
  3. 区分统计误差、数值与近似误差以及 Monte Carlo 误差；
  4. 根据计算目标，区分数值优化、数值积分和随机模拟的基本用途；
  5. 说明计算诊断、模型诊断和可复现报告为什么都是统计分析的必要环节；
  6. 从稳定性、准确性、计算成本和可复现性等方面评价一种计算方法。

本章示例给出可以按顺序直接运行的 R 代码；所有随机实验都显式设置种子，必要的扩展包也在代码中明确调用。读者暂时不必掌握全部语法，第 2 章将系统介绍 R 语言及可复现计算的基本方法。

1.2 统计建模的计算流程

一次完整的统计分析不只是“选择一个模型”或“调用一个函数”。研究问题、数据、推断目标、模型、算法和诊断共同决定了结论的含义与可信程度。从计算统计的角度看，统计建模可以概括为图 1.1 所示的流程。

统计建模的计算流程

图1.1: 统计建模的计算流程

  1. 研究问题与数据。明确所要回答的实质问题，检查数据来源、变量含义、缺失值、异常值和数据结构。
  2. 推断与预测目标。确定需要计算的对象，例如参数、标准误、置信区间、检验统计量、后验分布或预测分布。
  3. 模型。根据研究问题、目标和数据结构提出统计模型，例如线性模型、广义线性模型、层次模型或隐变量模型。
  4. 算法。选择可执行的计算方法，例如矩阵分解、数值优化、数值积分、Monte Carlo、Bootstrap 或 MCMC。
  5. 诊断。计算诊断检查算法是否收敛、数值结果是否稳定；模型诊断检查模型假设和拟合是否合理；预测诊断评价样本外表现是否可靠。
  6. 报告与复现。用可复现的代码、数据说明、图表和文字报告研究问题、模型、算法、诊断及结论。

1.2.1 以线性模型为例

考虑一个简单问题：解释变量 \(x\) 的变化是否与响应变量 \(y\) 的平均水平有关，并希望利用 \(x\) 预测新的\(y\)。若有 \(n\) 个观测 \((x_i,y_i)\)，可以从简单线性模型出发：

\[ y_i=\beta_0+\beta_1x_i+\varepsilon_i, \qquad i=1,\ldots,n. \]

这个例子虽然简单，却包含了统计建模流程中的各个环节，见表 1.1。

表1.1: 以简单线性模型说明统计建模的计算流程
环节 在线性模型中的具体内容 计算要点
研究问题与数据 研究 \(x\) 与 \(y\) 的平均关系；检查观测单位、缺失值、异常值和 \(x\) 的取值范围 数据定义决定系数和预测的解释范围
推断与预测目标 估计 \(\beta_0\)、\(\beta_1\)，量化斜率的不确定性，或预测给定 \(x\) 时的 \(y\) 参数估计、区间估计和预测是不同目标
统计模型 指定 \(E(y_i\mid x_i)=\beta_0+\beta_1x_i\)，并对误差的均值、方差及独立性作出假设 模型把实质问题转化为可估计的数学对象
原理与准则 普通最小二乘使残差平方和最小；在同方差正态误差下，它也给出 \(\beta\) 的极大似然估计 准则规定“什么结果算最好”
算法 构造设计矩阵 \(X\)，通常利用 QR 分解求解最小二乘问题 不必按公式显式计算 \((X^\top X)^{-1}\)
诊断 检查数值秩、条件数、残差结构、异常和高杠杆观测，预测任务还要检查样本外误差 算法稳定和模型合理需要分别验证
报告与复现 报告系数估计、标准误或区间、预测结果、诊断以及生成结果所需的代码和数据说明 数值输出必须结合假设和变量单位解释

因此，线性回归模型、最小二乘或似然原理、QR 分解和系数估计及预测结果处在不同层次。其中 \(\beta_0\)、\(\beta_1\) 是模型参数，利用样本计算出的 \(\hat\beta_0\)、\(\hat\beta_1\) 才是统计量。模型说明要估计什么，估计准则定义目标，算法把目标转化为有限步骤，诊断则判断这些步骤所得结果能否支持最终结论。后文将进一步说明，数学上可以写出的显式解，并不总是最适合计算机执行的形式。

图中的箭头并不表示这一过程只能单向进行。诊断可能促使研究者重新界定目标、修正模型或者更换算法；数据质量问题也可能要求重新整理数据甚至调整研究问题。统计分析因此是一个反复检查和修正的过程。

需要特别强调三个区别。第一，程序成功运行只说明计算过程没有被显式中断，并不保证结果数值可靠。第二，算法满足收敛条件并不意味着模型适合数据。第三，一项分析可以被完整复现，也仍然可能包含不恰当的假设或解释。可靠的统计结论需要统计理论、计算方法和实质判断共同支持。

1.3 从公式到计算

在课堂上推导统计公式时，我们通常把实数视为可以精确表示的对象，把代数运算视为可以无限精确地完成，并且把矩阵求逆、函数积分等符号看作可执行的数学运算。计算机并不在这个理想环境中工作。计算统计首先要求我们认识以下两个基本事实：

  1. 计算机中用于数值计算的浮点数只能表示实数的一个有限子集，浮点运算还会引入舍入误差，因而浮点算术并不完全遵循实数算术的代数性质 (Goldberg 1991; Higham 2002)；
  2. 数学上等价的表达式可能具有不同的数值性质，因此实际计算时需要根据舍入误差、问题的条件性和算法的稳定性选择合适的求值形式 (Wilkinson 1963; Higham 2002)。

计算机只能用有限位数表示数字。非常大的计算结果可能上溢，非常小的结果可能下溢，两个近似相等的数相减可能造成严重的有效数字损失。与此同时，直接翻译数学公式可能进行不必要的运算，也可能把一个原本可解的问题转化为数值上不稳定的问题。对于复杂似然函数、高维积分、隐变量模型和贝叶斯后验分布等没有简单解析解的对象，还必须依靠迭代或随机算法得到近似结果。这些现象及其误差分析是数值计算稳定性理论的基本内容 (Higham 2002)。

因此，面对一个具体的统计计算任务时，应当反复追问以下问题：

  1. 统计理论定义的目标量是什么？
  2. 哪一种算法能够稳定而有效地计算这个目标量？
  3. 数据规模扩大以后，时间和存储成本是否仍然可以接受？
  4. 统计误差、数值与近似误差以及模拟误差将如何影响结论？
  5. 应当使用什么诊断来判断结果是否可信？

1.4 例子：在对数尺度上计算

在统计建模中，似然函数往往是许多概率质量或密度贡献的乘积。即使每一项都可以正常表示，乘积也可能很快超出浮点数的表示范围。以标准正态分布为例：

dnorm(1)
#> [1] 0.242
dnorm(40)
#> [1] 0
dnorm(40, log = TRUE)
#> [1] -800.9

标准正态分布在 40 处的密度是一个极小的正数。dnorm(40) 返回 0，并不是因为打印时省略了小数位，而是因为结果在双精度浮点运算中发生了下溢。对数密度则仍然可以稳定表示。对于似然函数

\[ L(\theta)=\prod_{i=1}^n f(y_i \mid \theta), \]

实际计算通常使用对数似然

\[ \ell(\theta)=\log L(\theta)=\sum_{i=1}^n \log f(y_i \mid \theta). \]

对数变换避免了直接相乘大量极小数，并把乘法转化为加法，使求导、优化和程序实现更加方便。如果需要比较密度比

\[ \frac{f(x)}{g(x)}, \]

应先计算对数比

\[ \log f(x)-\log g(x). \]

当 \(f(x)>0\) 且 \(g(x)>0\) 时，只有在确实需要原尺度结果时，才计算 \(\exp\{\log f(x)-\log g(x)\}\)；如果目的只是比较两个量的大小，则应尽可能保留在对数尺度上，因为最后的指数运算仍然可能发生上溢或下溢。这个例子说明，数学上等价的表达式可能具有不同的数值表现 (Higham 2002; Monahan 2011)。

1.5 例子：线性回归不是直接求逆

线性回归是统计学中最常见的模型之一。令 \(y\in\mathbb R^n\) 为响应向量，\(X\in\mathbb R^{n\times p}\) 为设计矩阵，\(\beta\in\mathbb R^p\) 为未知参数。在经典同方差线性模型中，条件于 \(X\)，模型写作

\[ y=X\beta+\varepsilon. \]

并假设

\[ E(\varepsilon\mid X)=0, \qquad \operatorname{Var}(\varepsilon\mid X)=\sigma^2 I_n. \]

如果需要经典线性模型下的精确有限样本分布推断，通常进一步假设

\[ \varepsilon\mid X\sim N_n(0,\sigma^2 I_n). \]

正态性不是定义或计算普通最小二乘估计所必需的条件；它主要用于推导经典正态线性模型下的精确抽样分布、区间估计和检验。不同假设分别承担什么作用，是理解线性模型时需要特别区分的问题 (Seber and Lee 2003)。

普通最小二乘估计定义为

\[ \hat\beta \in \arg\min_{b\in\mathbb R^p}\lVert y-Xb\rVert_2^2. \]

任意最小二乘解都满足正规方程

\[ X^\top X\hat\beta=X^\top y. \]

如果 \(X\) 满列秩，即 \(\operatorname{rank}(X)=p\)，则最小二乘解唯一，\(X^\top X\) 可逆，因而可以在数学上写成

\[ \hat{\beta}=(X^\top X)^{-1}X^\top y. \]

这个公式适合揭示估计量的数学结构，但它并不意味着程序应当先求出 \((X^\top X)^{-1}\)。显式计算逆矩阵会完成求解 \(\hat\beta\) 所不需要的额外运算，并可能进一步放大舍入误差。较好的做法是直接求解正规方程。需要注意的是，这虽然避免了显式求逆，却仍然构造了 \(X^\top X\)。当 \(X\) 满列秩时，在二范数下有

\[ \kappa_2(X^\top X)=\kappa_2(X)^2, \]

因此正规方程会加重设计矩阵的病态程度。QR 分解不需要显式形成 \(X^\top X\)，通常是求解最小二乘问题更稳定的方法；在秩亏或近似秩亏的情形下，还可以使用秩揭示 QR 分解或奇异值分解。最小二乘问题的数值算法及稳定性分析可参见 Golub and Van Loan (2013) 和 Higham (2002)。

下面比较显式求逆、求解正规方程和 QR 分解。模拟中取 \(n=500\)、\(q=20\)，并令

\[ x_{ij}\stackrel{\mathrm{iid}}{\sim}N(0,1), \qquad i=1,\ldots,n,\quad j=1,\ldots,q, \]

其中 \(x_j=(x_{1j},\ldots,x_{nj})^\top\)，设计矩阵写作 \(X=(\boldsymbol 1_n,x_1,\ldots,x_q)\)。为生成一组固定的回归系数，先独立抽取

\[ \beta_0=1, \qquad \beta_j\stackrel{\mathrm{iid}}{\sim}\operatorname{Unif}(-1,1), \quad j=1,\ldots,q, \]

随后在给定这些系数的条件下生成

\[ y_i=\beta_0+\sum_{j=1}^{q}\beta_jx_{ij}+\varepsilon_i, \qquad \varepsilon_i\stackrel{\mathrm{iid}}{\sim}N(0,1). \]

这里随机抽取 \(\beta_1,\ldots,\beta_q\) 只是为了构造模拟数据；生成后将其视为固定参数。下表仅展示截距项以及与 \(x_1,\ldots,x_5\) 对应的计算结果。在设计矩阵条件良好时，三种方法给出的结果非常接近。

set.seed(2026)
n <- 500
p_slope <- 20
X_slope <- matrix(rnorm(n * p_slope), n, p_slope)
colnames(X_slope) <- paste0("x", seq_len(p_slope))
X <- cbind("(Intercept)" = 1, X_slope)
beta <- c(1, runif(p_slope, -1, 1))
sigma <- 1
y <- as.vector(X %*% beta + rnorm(n, sd = sigma))

ols_explicit_inverse <- function(X, y) {
  drop(solve(crossprod(X)) %*% crossprod(X, y))
}
ols_normal_equation <- function(X, y) {
  drop(solve(crossprod(X), crossprod(X, y)))
}
ols_qr <- function(X, y) {
  qr.coef(qr(X), y)
}

beta_inverse <- ols_explicit_inverse(X, y)
beta_normal <- ols_normal_equation(X, y)
beta_qr <- ols_qr(X, y)

coef_comparison <- data.frame(
  term = c("$1$", paste0("$x_", seq_len(5), "$")),
  explicit_inverse = beta_inverse[1:6],
  normal_equation = beta_normal[1:6],
  qr = beta_qr[1:6]
)
knitr::kable(
  coef_comparison,
  row.names = FALSE,
  digits = 5,
  col.names = c("模型项", "显式求逆", "正规方程", "QR 分解"),
  escape = FALSE
)
模型项 显式求逆 正规方程 QR 分解
\(1\) 1.02395 1.02395 1.02395
\(x_1\) -0.30539 -0.30539 -0.30539
\(x_2\) -0.33121 -0.33121 -0.33121
\(x_3\) 0.89234 0.89234 0.89234
\(x_4\) 0.37613 0.37613 0.37613
\(x_5\) 0.03115 0.03115 0.03115

计算时间也是算法比较的一部分。由于单次拟合很快，下面使用 microbenchmark 包在同一组 \(X\) 和 \(y\) 上交错重复运行三种方法，并报告耗时的中位数。相较于单次计时，多次交错运行可以减弱计时器分辨率、运行顺序和偶然系统负载的影响；使用中位数则可以降低少数异常耗时的影响。

library(microbenchmark)

timing_repetitions <- 300L
timing_benchmark <- suppressWarnings(
  microbenchmark(
    `显式求逆` = ols_explicit_inverse(X, y),
    `正规方程` = ols_normal_equation(X, y),
    `QR 分解` = ols_qr(X, y),
    times = timing_repetitions,
    unit = "ms",
    control = list(order = "random", warmup = 10L)
  )
)

timing_summary <- summary(timing_benchmark, unit = "ms")
timing_comparison <- data.frame(
  method = as.character(timing_summary$expr),
  median_ms = timing_summary$median,
  relative_time = timing_summary$median / min(timing_summary$median)
)
knitr::kable(
  timing_comparison,
  row.names = FALSE,
  digits = 3,
  col.names = c("方法", "耗时中位数（毫秒）", "相对耗时")
)
方法 耗时中位数（毫秒） 相对耗时
显式求逆 0.343 1.089
正规方程 0.315 1.000
QR 分解 0.550 1.745

具体耗时取决于硬件、线性代数库、系统负载和问题规模，因而不应把一次计时结果解释为普遍的速度排序。通常，显式求逆会完成求解系数并不需要的额外工作；直接求解正规方程往往更快，但会平方条件数；QR 分解兼顾稳定性，所需时间则可能高于正规方程。算法选择需要同时考虑计算时间、数值稳定性和问题结构，而不能只比较其中一项。

为了考察病态问题，在设计矩阵中加入一个与 \(x_1\) 几乎相同的变量。

set.seed(20260911L)
W <- cbind(X, x1_near_duplicate = X[, "x1"] + rnorm(n, sd = 1e-12))

qr_W <- qr(W)
normal_beta_W <- tryCatch(
  drop(solve(crossprod(W), crossprod(W, y))),
  error = function(e) rep(NA_real_, ncol(W))
)
qr_beta_W <- qr.coef(qr_W, y)

condition_summary <- data.frame(
  indicator = c("kappa(W)", "kappa(W^T W)", "QR 数值秩", "矩阵列数"),
  value = c(
    formatC(kappa(W), format = "e", digits = 3),
    formatC(kappa(crossprod(W)), format = "e", digits = 3),
    qr_W$rank,
    ncol(W)
  )
)
knitr::kable(
  condition_summary,
  row.names = FALSE,
  col.names = c("诊断量", "数值")
)
诊断量 数值
kappa(W) 3.344e+12
kappa(W^T W) 1.164e+16
QR 数值秩 21
矩阵列数 22 

ill_conditioned_coef <- data.frame(
  method = c("正规方程", "QR 分解"),
  x1 = c(normal_beta_W[2], qr_beta_W[2]),
  x1_near_duplicate = c(normal_beta_W[ncol(W)], qr_beta_W[ncol(W)])
)
knitr::kable(
  ill_conditioned_coef,
  row.names = FALSE,
  digits = 5,
  col.names = c("方法", "$x_1$", "$x_1$ 的近重复变量"),
  escape = FALSE
)
方法 \(x_1\) \(x_1\) 的近重复变量
正规方程 NA NA
QR 分解 -0.3054 NA

具体数值可能随计算平台和底层线性代数库略有差异。正规方程有时会直接报错，但更值得警惕的情形是：程序正常结束，却返回绝对值很大、方向相反并且对微小扰动十分敏感的系数。QR 分解把设计矩阵判断为数值秩不足，并用 NA 标出不能稳定区分的系数。这里的问题不仅是算法选择不当，也说明数据没有提供足够信息来分别识别两个近乎相同变量的作用。

这个例子体现了计算统计中的一个基本原则：程序没有报错并不等于结果可信。统计公式规定了目标，数值线性代数决定了实现路径，而条件数和数值秩等诊断帮助我们判断结果能否被可靠解释。需要注意，误差项分布假设决定的是模型的概率解释和统计推断；设计矩阵的条件数与数值秩决定的则是给定 \(X\) 和 \(y\) 后最小二乘问题能否被稳定计算。两类问题有关联，但不能相互替代。

1.6 例子：没有解析解时怎么办

统计公式可能把目标写得很清楚，却没有给出能够直接代入计算的解析表达式。不同目标需要不同的计算路径：求最优点通常使用数值优化，计算积分可以使用确定性的数值积分，也可以转化为随机模拟。本节用三个简短例子说明这些方法解决什么问题，以及计算结果应当怎样检查。

1.6.1 数值优化

在极大似然估计中，如果最优解位于参数空间内部并满足相应的正则条件，通常需要求解得分方程

\[ \frac{\partial \ell(\theta)}{\partial \theta}=0. \]

除少数简单模型外，这个方程往往不能手工求解。实际计算需要使用 Newton–Raphson、BFGS、梯度下降或随机梯度下降等算法迭代地寻找最优解。初值、停止准则、梯度和 Hessian 矩阵的计算精度都会影响最终结果；在非凸问题中，不同初值或算法还可能收敛到不同的局部最优解。因此，报告一个优化结果时，不能只给出参数估计，还应检查收敛状态、梯度大小、目标函数值以及对初值的敏感性。统计问题中的数值优化可参见 Monahan (2011) 和 Nocedal and Wright (2006)。

1.6.2 数值积分

另一类常见目标是计算定积分。例如，考虑

\[ I=\int_0^1 \exp(-x^2)\,dx. \]

这个积分没有初等函数形式的原函数，但这不妨碍我们计算它的数值近似。把区间 \([0,1]\) 等分为 \(m\) 段，令 \(h=1/m\)、\(x_j=jh\)，复合梯形公式给出

\[ \widehat I_m = h\left\{ \frac{q(x_0)+q(x_m)}{2} +\sum_{j=1}^{m-1}q(x_j) \right\}, \qquad q(x)=\exp(-x^2). \]

下面分别使用 100 段和 1000 段计算。这里的代码用于展示“网格加密”的思想，暂时不要求读者掌握函数定义的语法。

trapezoid_exp <- function(m) {
  x_grid <- seq(0, 1, length.out = m + 1)
  q_grid <- exp(-x_grid^2)
  h <- 1 / m
  h * (sum(q_grid) - (q_grid[1] + q_grid[m + 1]) / 2)
}

integral_100 <- trapezoid_exp(100)
integral_1000 <- trapezoid_exp(1000)

c(
  grid_100 = integral_100,
  grid_1000 = integral_1000,
  refinement_difference = abs(integral_1000 - integral_100)
)
#>              grid_100             grid_1000 refinement_difference 
#>             7.468e-01             7.468e-01             6.070e-06

两次计算使用确定的网格，因此在相同输入和计算环境下会得到相同结果。把网格从 100 段加密到 1000 段后，两次近似已经十分接近；这种比较可以作为误差检查，但两次结果接近本身并不能证明它们等于精确值。数值积分通常还会根据误差估计自动调整网格，用户则通过绝对或相对容差控制所需精度。第 7 章将系统介绍这些方法。

1.6.3 随机模拟

有些积分可以写成期望，并通过随机抽样近似。若目标量为 \(\mu=E\{g(X)\}\)，从 \(X\) 的分布中生成\(M\) 个样本 \(X_1,\ldots,X_M\) 后，可以用样本平均值

\[ \widehat\mu_M= \frac{1}{M}\sum_{i=1}^M g(X_i), \]

近似 \(\mu\)。这就是 Monte Carlo 方法的基本思想。前面的积分可以写成

\[ I=E\{\exp(-U^2)\}, \qquad U\sim \operatorname{Unif}(0,1). \]

因此，只需从 \([0,1]\) 均匀抽样，再计算 \(\exp(-U^2)\) 的平均值：

set.seed(123)
M <- 1e4
u_mc <- runif(M)
g_mc <- exp(-u_mc^2)

c(
  estimate = mean(g_mc),
  estimated_mcse = sd(g_mc) / sqrt(M),
  trapezoid_reference = integral_1000
)
#>            estimate      estimated_mcse trapezoid_reference 
#>            0.748874            0.001994            0.746824

Monte Carlo 结果会随随机样本改变，因此代码应设置随机种子，并用 Monte Carlo 标准误（MCSE）描述模拟误差。上面用样本标准差除以 \(\sqrt M\) 估计 MCSE；模拟次数越多，结果通常越稳定。对于这个光滑的一维积分，梯形公式通常更有效；随机模拟的价值主要体现在高维或结构复杂的问题中。第 10 章将系统介绍 Monte Carlo 方法。

1.6.4 计算目标决定计算方法

前面的例子说明，“没有解析解”只是计算问题的起点。首先要识别数学目标，再选择相应的算法和诊断：

计算目标 常用计算方法 主要输出 常用检查方法
解线性方程或最小二乘问题 矩阵分解 解或系数估计 残差、条件数、数值秩
求方程的零点 数值求根 近似根 方程残差、区间或初值敏感性
求函数的最大值或最小值 数值优化 近似最优点 收敛码、梯度、目标函数值、初值敏感性
计算低维光滑积分 数值积分 确定性近似值 误差估计、容差、网格加密
计算高维积分或复杂期望 Monte Carlo 随机近似值 模拟次数、Monte Carlo 标准误

表中的对应关系不是固定配方。同一个积分既可以用求积公式，也可以用 Monte Carlo；得分方程既可以直接求根，也可以通过优化对数似然求解。方法的选择还取决于维数、函数是否光滑、能否计算导数、可用计算资源以及所需精度。同一个统计分析也可能连续用到多种方法：先通过矩阵分解求解线性方程，再用数值优化估计参数，最后用随机模拟评价某个难以直接积分的期望。

本书后续各部分将分别系统介绍这些工具，并把它们用于回归、重采样及其他统计问题。第五部分学习贝叶斯推断时，读者还会看到数值优化、数值积分和随机模拟如何围绕一种新的推断目标组合起来；相关概率公式将在那里正式引入，此处不要求预先掌握。上面的网格误差和 Monte Carlo 误差也提示我们，不同近似方法会带来不同类型的误差，下一节将对此作进一步区分。

1.7 计算统计中的三类误差

计算结果与目标值之间存在差异，并不必然意味着模型错误或程序存在缺陷。为了正确解释计算结果，需要区分误差的来源。

第一类是统计误差。即使计算完全精确，有限观测样本也不能提供总体的全部信息。样本均值与总体均值、参数估计与真实参数之间的差异来自数据的随机性。标准误、置信区间和后验不确定性主要刻画这类误差。增加有效的观测信息通常可以降低统计误差。

第二类是数值与确定性近似误差。浮点数舍入、上溢、下溢和病态矩阵会产生数值误差；迭代算法过早停止、积分网格过粗或级数截断则会产生确定性近似误差。条件数、线性方程残差、目标函数变化、梯度大小以及不同容差下的敏感性分析，都是评价这类误差的重要工具。简单地增加观测样本量并不一定能够改善数值稳定性，有时反而会使计算规模和病态程度增加。

第三类是模拟误差。使用 Monte Carlo、Bootstrap 或 MCMC 方法时，有限模拟次数会使结果随随机数种子而变化。独立 Monte Carlo 中的误差通常按 \(M^{-1/2}\) 的速度下降；MCMC 样本存在自相关，因而还需要考虑有效样本量。模拟结果不仅应报告点估计，也应报告 Monte Carlo 标准误或其他能够反映计算精度的指标 (Robert and Casella 2004)。

这三类误差对应不同的控制方法。增加观测样本量主要影响统计误差，改变数值算法和容差主要影响数值与近似误差，增加有效模拟次数主要影响模拟误差。混淆这些误差可能导致无效的补救措施。例如，增加 MCMC 迭代次数不能弥补模型设定错误，也不能消除实际数据样本量不足造成的统计不确定性。

模型设定偏差不属于上述三类计算误差，但它同样可能使统计结论失真。残差分析、预测检验、样本外评价和敏感性分析用于考察模型与数据及研究问题是否相符。算法收敛是得到可信结果的必要条件之一，却不是模型正确的充分条件。

1.8 计算统计的基本问题

从以上例子可以看出，计算统计至少涉及以下五类基本问题。

第一，数值表示与稳定性。计算机中的浮点数只能有限精度地表示实数。浮点误差、上溢、下溢、灾难性抵消和病态矩阵都可能改变计算结果。对数尺度、标准化、矩阵分解和条件数分析是处理这些问题的基本工具。

第二，确定性算法。求根、优化、数值积分、矩阵分解和 EM 算法等方法，使没有简单解析解的估计与近似问题成为可计算的问题。算法不仅要产生结果，还应提供误差控制或收敛诊断。

第三，随机模拟。当积分、分布或不确定性无法直接计算时，可以借助随机数生成、Monte Carlo、重采样和 MCMC 方法构造近似推断，并通过 Monte Carlo 标准误或有效样本量评价模拟精度。

第四，效率与可扩展性。矩阵规模、参数维数和模拟次数决定了时间与存储成本。一个在小数据上可行的方法，在大样本或高维问题中可能无法运行。计算复杂度、稀疏结构、向量化和并行计算都会影响算法选择。

第五，验证与可复现性。统计分析不仅要算得出，还要能够被检查、复现和扩展。完整的代码、数据来源、随机种子、软件及软件包版本、诊断结果和分析说明，都是现代统计实践的重要组成部分。

1.9 各部分的关系

全书的分部没有逐项对应上一节的五类问题。章节的先后主要取决于后面的内容要用到什么：先讲 R 和数值计算基础，再分别学习确定性算法和随机模拟，然后把这些方法用于统计模型，再进入贝叶斯计算。图 1.2 画出前五部分的关系，最后的现代应用部分综合使用这些方法。

本书前五部分的关系。实线表示正文顺序，虚线表示贝叶斯计算会直接用到确定性方法和随机模拟。

图1.2: 本书前五部分的关系。实线表示正文顺序，虚线表示贝叶斯计算会直接用到确定性方法和随机模拟。

第一部分先处理计算的基本问题。本章说明一个统计目标怎样转化成计算任务；第 2 章介绍 R 语言和可复现计算；第 3 章讨论浮点数、误差传播和数值稳定性。把这些内容放在前面，是因为后面几乎每一种方法都要和误差打交道。如果分不清统计误差、数值与近似误差和模拟误差，即使程序顺利运行，也很难判断结果是否可信。

第二部分讨论在没有简单解析解时怎样做确定性计算。矩阵方法排在最前，是因为后面的优化、Laplace 近似和 EM 计算都会反复用到线性方程与矩阵分解；其后依次讨论求根、优化、数值积分和 EM 算法。输入、初值和计算设置不变时，确定性算法不会因随机抽样而产生波动；但仍需检查数值稳定性、近似误差和迭代收敛情况。

第三部分转向随机数、Monte Carlo 和重采样。Monte Carlo 用模拟样本近似积分和概率，Bootstrap、Jackknife 和置换检验则从已有数据出发考察统计量的抽样变化。使用随机模拟时，除了数据本身带来的统计误差，还要说明有限模拟次数给结果增加了多大误差。

到了第四部分，重点从“一个算法怎样工作”转到“一个模型怎样算”。最小二乘依赖矩阵分解，广义线性模型需要迭代求解，交叉验证要反复划分数据，正则化则离不开优化。把它们放在基本算法之后，是为了把计算步骤和统计解释放在一起看：算法停下来了，不等于模型合理；样本内拟合得好，也不等于样本外预测可靠。

第五部分进入贝叶斯计算。写出先验、似然和后验只是确定了目标；后验众数常要通过优化求得，没有解析解的后验积分要用数值积分、Laplace 近似或 Monte Carlo，回归模型还要用到矩阵方法。这里并不是另起炉灶，而是让前面的工具围绕后验分布重新组合。

第六部分是第 20 章的综合应用。生成模型中的噪声回归、文本采样的概率修正、水印检测与预测驱动推断，分别承接优化、随机模拟和统计推断；低秩适配与偏好优化则回到矩阵计算和二元似然。读者可以通过这些案例检验自己能否把已有方法用于新的计算问题。

效率与可扩展性、验证与可复现性没有单独占一个部分，因为每一章都要面对它们：数据或参数规模变大以后还能不能算，得到的结果能不能检查和复现。第一次阅读沿图中的实线即可；已有 R 和数值计算基础的读者可以快读第一部分，但最好不要略过本章关于三类误差的区分。

MCMC 看起来本来应该放在第三部分：从算法分类看，它确实是一种随机模拟方法。本书仍把它放在贝叶斯部分，主要因为这里用它解决后验抽样问题。如果还没有见过后验计算的目标，也不知道网格法、Laplace 近似和重要性抽样会在哪里遇到困难，MCMC 很容易只剩下一套抽样步骤。因此，读者先学习第 15 章和第 16 章，再进入第 17 章。这只是讲授次序的选择，MCMC 并不限于贝叶斯问题。

1.10 进一步阅读

关于计算统计的整体框架，可阅读 Gentle (2009) 和 Monahan (2011)；关于浮点运算、误差传播与数值稳定性，可阅读 Higham (2002)；关于矩阵分解与最小二乘计算，可阅读 Golub and Van Loan (2013)；关于数值优化，可阅读 Nocedal and Wright (2006)；关于 Monte Carlo 方法及其误差分析，可阅读 Robert and Casella (2004)；关于贝叶斯建模、后验计算和模型检查，可阅读 Gelman et al. (2013)。线性回归模型的概率假设、估计与推断可参见 Seber and Lee (2003)。

1.11 本章小结

本章的核心观点是：统计理论定义推断目标，计算统计提供可执行路径。数学上正确的表达式未必适合直接交给计算机求值，能够运行的程序也未必产生稳定可信的结果。可靠的计算实现需要同时关注有限精度、算法选择、计算成本、误差评估、模型诊断和可复现性。

对数尺度计算说明了数学等价不代表数值等价；病态回归说明了显式求逆、正规方程和 QR 分解在计算时间与数值稳定性上各有差异，也说明“无报错”不能替代数值诊断；Monte Carlo 近似则说明随机模拟在产生估计的同时还会引入可以量化的计算误差。没有解析解时，需要先辨认目标是求最优点、算积分还是生成样本，再选择数值优化、数值积分或随机模拟方法。后续各章将进一步发展这些思想，使读者不仅知道一个统计量的定义，还能够判断它应当如何计算以及计算结果是否可信。

1.12 思考题

  1. 结合本章的例子说明，为什么数学上等价的表达式未必具有相同的数值表现？“程序没有报错”为什么不能作为结果可信的充分证据？
  2. 在线性回归中，分别说明下列条件的作用：\(X\) 满列秩、\(E(\varepsilon\mid X)=0\)、\(\operatorname{Var}(\varepsilon\mid X)=\sigma^2I_n\) 以及 \(\varepsilon\mid X\) 服从正态分布。哪些条件关系到最小二乘解的唯一性，哪些条件关系到估计量性质或经典有限样本推断？
  3. 比较显式求逆、求解正规方程和 QR 分解三种最小二乘计算方法。为什么 solve(crossprod(X), crossprod(X, y)) 虽然优于显式求逆，却仍可能比 QR 分解更容易受到病态矩阵影响？
  4. 对下列情形分别判断主要涉及统计误差、数值与近似误差、模拟误差还是模型设定偏差，并说明应当采用什么诊断或改进方法：观测样本量过小；设计矩阵接近奇异；MCMC 有效样本量过低；逻辑回归遗漏重要的非线性关系；数值积分网格过粗。
  5. 分别解释“算法已经收敛”“模型能够合理描述数据”和“分析可以被复现”的含义。为什么三者不能相互替代？
  6. 对于同一个一维积分，数值积分和 Monte Carlo 都可以给出近似。比较两种方法的结果是否随随机数种子变化、怎样提高精度以及应报告什么误差或诊断；再说明在什么情况下可能更适合使用 Monte Carlo。
  7. 选择一个熟悉的经济或社会统计问题，依次说明研究问题与数据、推断或预测目标、统计模型、计算算法、计算诊断、模型诊断和可复现报告应当包含哪些内容。

1.13 上机实验（Lab）

上机实验的目标不仅是得到一个数值结果，还要形成可以复核的计算证据。每份实验报告至少应包括问题与目标、可运行代码、随机种子、主要输出或图表、必要的计算诊断以及对结果的解释。

  1. Lab 1：对数尺度与浮点下溢。 对标准正态分布，比较 dnorm(10)、dnorm(40) 和 dnorm(40, log = TRUE) 的输出，解释下溢与输出位数截断的区别。再设两个对数密度为 \(-1000\) 和 \(-1001\)，分别在原尺度和对数尺度上计算或比较密度比，并说明何时没有必要执行指数运算。
  2. Lab 2：病态最小二乘问题。 在本章的回归示例中，把近重复变量的噪声标准差依次改为 \(10^{-4}\)、\(10^{-8}\) 和 \(10^{-12}\)。比较 \(\kappa_2(X)\)、\(\kappa_2(X^\top X)\)、QR 数值秩和两个相关变量的系数。再对响应变量施加一个很小的扰动，分别考察回归系数、拟合值和残差的变化，并解释为什么三者的敏感程度可能不同。
  3. Lab 3：确定性积分与随机积分。 对 \(I=\int_0^1\exp(-x^2)\,dx\)，先把梯形公式的分段数依次设为 \(m=10,100,1000,10000\)，以 R 的 integrate() 结果为参照考察误差怎样变化；再令 \(M=100,1000,10000,100000\)，各重复 500 次 Monte Carlo 估计，计算经验标准差和平均 MCSE。比较两类近似的结果是否随随机数种子变化、怎样提高精度以及各自需要记录的计算设置。
  1. 拓展 Lab：多元正态对数密度。 利用 Cholesky 分解计算 \(p\) 维正态对数密度中的二次型与对数行列式，不显式求协方差矩阵的逆或行列式。用模拟数据比较该实现与直接计算在数值结果和用时上的差异，并说明其适用条件。

参考文献

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. Boca Raton, FL: CRC Press.
Gentle, James E. 2009. Computational Statistics. New York: Springer.
Goldberg, David. 1991. “What Every Computer Scientist Should Know about Floating-Point Arithmetic.” ACM Computing Surveys 23 (1): 5–48. https://doi.org/10.1145/103162.103163.
Golub, Gene H., and Charles F. Van Loan. 2013. Matrix Computations. 4th ed. Baltimore: Johns Hopkins University Press.
Higham, Nicholas J. 2002. Accuracy and Stability of Numerical Algorithms. 2nd ed. Philadelphia: SIAM.
Monahan, John F. 2011. Numerical Methods of Statistics. 2nd ed. Cambridge: Cambridge University Press.
Nocedal, Jorge, and Stephen J. Wright. 2006. Numerical Optimization. 2nd ed. New York: Springer.
Peng, Roger D. 2018. “Advanced Statistical Computing.” Work in Progress. https://bookdown.org/rdpeng/advstatcomp/.
Robert, Christian P., and George Casella. 2004. Monte Carlo Statistical Methods. 2nd ed. New York: Springer.
Seber, George A. F., and Alan J. Lee. 2003. Linear Regression Analysis. 2nd ed. Hoboken, NJ: Wiley.
Wilkinson, J. H. 1963. Rounding Errors in Algebraic Processes. London: Her Majesty’s Stationery Office.