第 2 章 R 语言与可复现计算
第 1 章已经说明，计算统计关心的是如何把统计目标转化为可靠、稳定、高效且可验证的计算过程。本章介绍本书后续章节所需的 R 语言基础，并把重点放在可复现计算上。这里的目标不是完整讲授 R 语言，而是建立一种适合计算统计的工作方式：用清晰的数据对象表示问题，用函数封装计算步骤，用随机数种子控制模拟实验，用项目目录和动态文档保留完整的计算证据。
本章中的代码只使用少量常用函数。没有 R 编程经验的读者可以按顺序运行示例；已有经验的读者也应关注本章关于项目组织、随机数控制、函数化计算和报告复现的原则，因为这些原则会贯穿后续的数值线性代数、优化、Monte Carlo、重采样和贝叶斯计算。
2.1 本章学习目标
完成本章后，读者应当能够：
- 说明 R 在计算统计中的作用，并掌握基本的帮助、包安装和项目组织方式；
- 使用向量、矩阵、列表和数据框表示常见统计数据；
- 正确使用子集提取、向量化和矩阵运算，避免常见的维度错误；
- 用控制流和函数封装可重复执行的计算任务；
- 使用
set.seed()、相对路径、R Markdown 和环境记录提高分析的可复现性； - 根据变量类型选择合适的基础图形或
ggplot2图形，并理解图形在统计计算诊断中的作用。
2.2 R 在计算统计中的作用
R 是一种面向统计计算和统计图形的语言与环境。它继承了 S 语言的许多思想，并由 Ross Ihaka 和 Robert Gentleman 在 20 世纪 90 年代发展为开源系统；关于 R 的早期设计和 S 语言背景，可参见 Ihaka and Gentleman (1996) 和 Chambers (2008)。R 的优势不只在于语法本身，还在于围绕统计建模、数值算法、图形展示、动态文档和软件包形成的生态系统。
本书采用 R 作为主要计算语言，原因有三点。第一，许多统计方法会先以 R 包的形式发布，R 因此是学习和复现实证方法的重要平台。第二，R 的向量、矩阵和数据框与统计问题中的基本对象高度一致，适合直接表达模型矩阵、样本、参数和模拟结果。第三，R Markdown、knitr、bookdown、Quarto 等工具使代码、公式、图表、文字解释和参考文献可以放在同一份动态文档中，这对计算统计中的复现和审稿尤其重要；动态文档和可复现研究的思想可参见 Gentleman and Temple Lang (2007)、Roger D. Peng (2011)、Xie (2015) 和 Xie, Allaire, and Grolemund (2018)。
本章将 R 视为“计算统计的工作台”。后续章节讨论的线性方程求解、矩阵分解、数值优化、随机模拟和 MCMC 方法，都会通过 R 代码展示其计算行为。
2.3 工作环境与项目组织
R 可以从 CRAN 网站下载安装；RStudio IDE 或 Positron 等集成开发环境可以提供代码编辑、对象查看、图形预览、调试和文档渲染等功能。安装外部包时，可以使用
安装完成后，用 library() 将包载入当前 R 会话。例如，本章后面使用 ggplot2 作图时会运行 library(ggplot2)。
帮助系统
R 的帮助系统是学习函数最可靠的入口。若要查看函数 solve() 的帮助文档，可以运行
或
如果只记得关键词而不确定函数名，可以使用
函数帮助页通常包含函数用途、参数含义、返回值、注意事项和示例。学习一个新函数时，不应只看示例代码，还应阅读参数说明，尤其是默认值、缺失值处理和返回对象的结构。
项目目录
可复现分析应当以项目为单位组织。一个项目通常包含数据、代码、输出、报告和说明文件，例如：
project/
data-raw/ 原始数据
data/ 清理后的分析数据
R/ 函数和脚本
output/ 表格、图形和中间结果
report.Rmd 动态报告
README.md 项目说明项目目录的核心原则是：代码中使用相对于项目根目录的路径，而不是依赖某台电脑上的绝对路径。这样，项目移动到另一台机器或交给其他人时，仍然可以运行。
在 R 中，可以用 getwd() 查看当前工作目录：
交互式分析中可以临时使用 setwd() 改变工作目录，但在教材、论文复现代码和正式项目中，不建议把机器相关的绝对路径写入脚本。更稳妥的方式是使用 RStudio Project、Quarto Project 或其他项目根目录机制，并保持数据与代码的相对路径关系。
2.4 对象与基本数据结构
R 中的一切计算都围绕对象展开。一个对象可以是数值、向量、矩阵、列表、数据框、函数或模型拟合结果。使用 str() 可以查看对象的内部结构，这在理解函数返回值时特别重要。
向量
向量是 R 中最基本的数据结构。同一个原子向量中的元素必须属于同一种基本类型，常见类型包括逻辑型、整数型、数值型和字符型。
num_x <- c(1.2, 4.5, 6)
int_x <- c(1L, 6L, 80L)
log_x <- c(TRUE, FALSE, TRUE)
cha_x <- c("mean", "variance", "likelihood")
str(num_x)
#> num [1:3] 1.2 4.5 6
str(int_x)
#> int [1:3] 1 6 80
str(log_x)
#> logi [1:3] TRUE FALSE TRUE
str(cha_x)
#> chr [1:3] "mean" "variance" "likelihood"R 提供了若干生成规则向量的函数。
10:15
#> [1] 10 11 12 13 14 15
seq(from = 1, to = 10, by = 3)
#> [1] 1 4 7 10
rep(14, times = 5)
#> [1] 14 14 14 14 14
LETTERS[1:5]
#> [1] "A" "B" "C" "D" "E"统计数据中经常出现缺失值。R 使用 NA 表示缺失。许多函数在遇到 NA 时会返回 NA，除非明确要求忽略缺失值。
缺失值处理是统计分析的一部分。简单地使用 na.rm = TRUE 可以让程序继续运行，但并不自动解决缺失机制带来的统计问题。
矩阵
矩阵是计算统计中最常见的对象之一。模型矩阵、协方差矩阵、Hessian 矩阵和转移矩阵都可以用矩阵表示。
X <- matrix(1:6, nrow = 2, ncol = 3)
X
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6
dim(X)
#> [1] 2 3
t(X)
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 3 4
#> [3,] 5 6行合并和列合并分别使用 rbind() 与 cbind()。
数据框
数据框是一种二维数据结构，每一列可以有不同类型，但各列长度必须相同。它是 R 中表示观测数据的主要结构。
dat <- data.frame(
id = 1:5,
group = c("A", "A", "B", "B", "B"),
y = c(2.1, 2.4, 3.0, 2.8, 3.2)
)
str(dat)
#> 'data.frame': 5 obs. of 3 variables:
#> $ id : int 1 2 3 4 5
#> $ group: chr "A" "A" "B" "B" ...
#> $ y : num 2.1 2.4 3 2.8 3.2
head(dat)
#> id group y
#> 1 1 A 2.1
#> 2 2 A 2.4
#> 3 3 B 3.0
#> 4 4 B 2.8
#> 5 5 B 3.2数据框可以看作特殊的列表：每一列是列表中的一个元素。正因为如此，数据框既支持矩阵式的行列提取，也支持列表式的列提取。
2.5 子集提取、向量化与矩阵运算
R 中常用的子集提取操作符有三种：
[返回与原对象相同类型的子集；[[从列表或数据框中提取单个元素；$用名称从列表或数据框中提取元素。
下面的例子展示三者的差异。
result <- list(
coefficient = c(1.2, -0.7),
vcov = matrix(c(0.4, 0.1, 0.1, 0.3), 2, 2),
converged = TRUE
)
result[1]
#> $coefficient
#> [1] 1.2 -0.7
result[[1]]
#> [1] 1.2 -0.7
result$coefficient
#> [1] 1.2 -0.7对矩阵和数据框提取子集时，要特别注意维度是否被自动下降。下面的第一行返回向量，第二行保留矩阵结构。
X <- matrix(1:9, nrow = 3)
X[, 1]
#> [1] 1 2 3
X[, 1, drop = FALSE]
#> [,1]
#> [1,] 1
#> [2,] 2
#> [3,] 3在计算统计中，维度错误常常比语法错误更隐蔽。例如，一个函数理论上要求输入矩阵，但某次子集提取返回了向量，程序可能仍然运行，却改变了后续计算的含义。因此，编写函数时应在关键位置检查输入维度。
R 的许多运算是向量化的。向量化不仅能使代码更简洁，也通常能提高计算效率。
矩阵运算应尽量使用专门函数。例如，crossprod(X) 等价于 t(X) %*% X，但通常更清晰，也可能更高效。
X <- cbind(1, 1:4, c(2, 1, 3, 5))
beta <- c(0.5, 1.0, -0.2)
X %*% beta
#> [,1]
#> [1,] 1.1
#> [2,] 2.3
#> [3,] 2.9
#> [4,] 3.5
crossprod(X)
#> [,1] [,2] [,3]
#> [1,] 4 10 11
#> [2,] 10 30 33
#> [3,] 11 33 39需要注意，向量化并不意味着所有循环都应避免。更重要的是让代码表达清楚，并避免在循环中反复扩展对象长度。若必须循环，通常应预先分配结果对象。
2.6 控制流、循环与函数
控制流用于根据条件决定执行路径，循环用于重复执行某段代码，函数用于把一段计算封装成可复用的单位。
循环函数
R 中也常用 lapply()、sapply()、vapply()、apply() 和 mapply() 表达重复计算。与 sapply() 相比，vapply() 要求事先指定返回值类型，因此在正式函数中更稳妥。
set.seed(123)
samples <- list(
normal = rnorm(100),
uniform = runif(100),
t_3 = rt(100, df = 3)
)
vapply(samples, mean, numeric(1))
#> normal uniform t_3
#> 0.09041 0.48605 0.02283apply() 可用于矩阵的行或列。
编写函数
当一段代码需要重复运行、需要更换参数、或者代表一个明确的计算步骤时，就应当写成函数。函数的基本形式为：
下面的函数计算正态样本在参数向量\(\boldsymbol\theta=(\mu,\log\sigma)^\top\) 下的对数似然。这里使用 \(\log\sigma\) 而不是直接使用\(\sigma\)，是为了保证标准差为正。
normal_loglik <- function(theta, y) {
if (length(theta) != 2L) {
stop("theta must contain mu and log_sigma.")
}
mu <- theta[1]
sigma <- exp(theta[2])
sum(dnorm(y, mean = mu, sd = sigma, log = TRUE))
}
set.seed(2026)
y <- rnorm(50, mean = 2, sd = 1.5)
theta0 <- c(mu = mean(y), log_sigma = log(sd(y)))
normal_loglik(theta0, y)
#> [1] -89.5这个例子体现了计算统计中编写函数的几个基本习惯：输入参数有明确含义；非法输入应尽早报错；概率密度计算优先使用对数尺度；函数返回一个可供后续优化或比较使用的数值。
2.7 随机数与可复现模拟
随机模拟是计算统计的重要组成部分。R 中许多分布函数具有统一的命名规则。以正态分布为例：
| 前缀 | 含义 | 示例 |
|---|---|---|
| d | 密度或概率质量函数 | dnorm() |
| p | 分布函数 | pnorm() |
| q | 分位数函数 | qnorm() |
| r | 随机数生成 | rnorm() |
同一段随机模拟代码，如果没有控制随机数种子，每次运行一般会得到不同结果。使用 set.seed() 可以把随机数生成器置于一个确定状态，从而复现实验结果。
set.seed(2026)
rnorm(5)
#> [1] 0.52059 -1.07969 0.13924 -0.08475 -0.66664
set.seed(2026)
rnorm(5)
#> [1] 0.52059 -1.07969 0.13924 -0.08475 -0.66664在模拟研究中，通常应在实验开始前设置一次种子，然后让随机数序列自然向前推进。下面的函数估计\(\operatorname{E}(X^2)\)，其中 \(X\sim N(0,1)\)。理论值为 1。
estimate_ex2 <- function(n, seed = NULL) {
if (!is.null(seed)) {
set.seed(seed)
}
mean(rnorm(n)^2)
}
estimate_ex2(10000, seed = 1)
#> [1] 1.025
estimate_ex2(10000, seed = 1)
#> [1] 1.025如果在重复实验的每一次都使用同一个种子，就会得到完全相同的重复结果，这通常不是我们想要的 Monte Carlo 重复实验。
replicate(3, estimate_ex2(1000, seed = 1))
#> [1] 1.07 1.07 1.07
set.seed(1)
replicate(3, estimate_ex2(1000))
#> [1] 1.070 1.081 1.062第一行代码重复了同一个实验三次；第二种写法则生成三个独立的模拟重复。随机数种子不是形式细节，它直接决定模拟实验是否可复现、是否真正重复。
2.8 数据输入输出
计算统计中的数据输入输出应服务于可复现流程。原始数据应尽量保留不变，清洗过程应写成脚本或函数，中间结果应有明确来源。
常用的数据输入函数包括：
read.csv()、read.table()：读取表格数据；readLines()：读取文本行；readRDS()：读取单个 R 对象；source()：运行 R 脚本。
常用的数据输出函数包括：
write.csv()、write.table()：输出表格数据；writeLines()：输出文本行；saveRDS()：保存单个 R 对象；sink()：将控制台输出临时写入文件。
下面用临时文件演示表格数据的读写。
sim_data <- data.frame(
x = 1:5,
y = c(2.0, 2.3, 2.7, 3.1, 3.6)
)
tmp_csv <- tempfile(fileext = ".csv")
write.csv(sim_data, tmp_csv, row.names = FALSE)
read.csv(tmp_csv)
#> x y
#> 1 1 2.0
#> 2 2 2.3
#> 3 3 2.7
#> 4 4 3.1
#> 5 5 3.6若需要保存一个 R 对象，saveRDS() 和 readRDS() 通常比 save() 和 load() 更适合可复现流程，因为它们读写的是单个对象，赋值名称由代码明确给出。
tmp_rds <- tempfile(fileext = ".rds")
saveRDS(sim_data, tmp_rds)
sim_data_again <- readRDS(tmp_rds)
identical(sim_data, sim_data_again)
#> [1] TRUE相比之下，load() 会把文件中保存的对象名直接恢复到当前环境中。如果文件中对象较多，或者对象名称与当前环境已有对象冲突，可能降低代码的可读性和可控性。
2.9 可复现报告与计算环境
可复现计算要求读者能够从同一份数据和代码重新得到相同或可解释差异范围内的结果。对于计算统计而言，这一要求至少包括五类信息：
- 数据来源、清理规则和样本筛选条件；
- 执行计算的代码和参数设置；
- 随机数种子和随机数生成方法；
- R 版本、包版本和操作系统等环境信息；
- 图表、表格和结论之间的对应关系。
R Markdown 是组织这些信息的常用工具。一个最小的 R Markdown 报告可以包含标题、正文和代码块。下面的示意中，[r] 表示一个 R 代码块；在实际 R Markdown 文件中，应将其写为三个反引号后接 {r}。
---
title: "计算统计实验报告"
output: html_document
---
本报告用一个简单数据集演示描述统计。
```[r]
summary(cars)
plot(cars)
```渲染报告时，代码会被重新执行，输出、图形和文字会合并成 HTML、PDF 或 Word 等格式。这样可以减少“图是手工复制的、表是旧版本结果”的风险。
正式分析中，建议在报告末尾记录环境信息：
对于需要长期维护或多人协作的项目，可以使用 renv 等工具记录包版本。即便不使用专门工具，也应在报告或 README 中说明主要软件版本、数据版本和生成结果的命令。
2.10 数据可视化
图形不是对结果的装饰，而是统计计算的重要诊断工具。模型拟合前，图形帮助我们理解数据结构、异常值和变量关系；模型拟合后，图形帮助我们检查残差、收敛过程、模拟误差和预测表现。
变量类型决定了图形类型。条形图和饼图通常用于分类变量；直方图、密度图和箱线图通常用于定量变量；散点图用于展示两个定量变量之间的关系。R 的 graphics 包提供了基础绘图函数，常见函数如表 2.2 所示。
| R.函数 | 主要用途 |
|---|---|
barplot() |
分类变量的数量或比例 |
pie() |
分类变量的组成比例 |
hist() |
定量变量的频数分布 |
density() |
定量变量的平滑分布 |
boxplot() |
定量变量的分位数和异常值 |
plot() |
通用绘图函数 |
smoothScatter() |
高密度散点图 |
pairs() |
多变量两两散点图 |
image() |
矩阵或网格数据图像 |
contour() |
等高线图 |
persp() |
三维曲面图 |
基础图形
以居民消费支出为例，分类变量“类别”对应不同消费项目，条形图可以比较不同类别的支出水平。
expense <- data.frame(
spending = c(6084, 5055, 2862, 2513, 1902, 1338, 1281, 524),
category = c(
"食品烟酒", "居住", "交通通信", "教育文化娱乐",
"医疗保健", "衣着", "生活用品及服务", "其他用品及服务"
)
)
par(family = cn_plot_family)
barplot(
height = expense$spending,
names.arg = expense$category,
main = "居民人均消费支出示例",
ylab = "消费（元）",
border = "darkblue",
col = "orange3",
las = 2
)
直方图用于展示定量变量的分布。下面生成一个课程成绩示例。
set.seed(123)
grades <- sample(70:100, 50, replace = TRUE)
par(family = cn_plot_family)
hist(
grades,
breaks = seq(70, 100, by = 5),
col = "orange3",
border = "darkblue",
main = "《计算统计》课程成绩直方图",
xlab = "成绩",
ylab = "频数"
)
箱线图基于分位数展示定量变量的分布。设第一四分位数为 \(Q_1\)，第三四分位数为 \(Q_3\)，四分位距为 \(\mathrm{IQR}=Q_3-Q_1\)。通常把 \(Q_1-1.5\mathrm{IQR}\) 和 \(Q_3+1.5\mathrm{IQR}\) 称为下、上界限；箱线图的须端画到落在界限内的最小和最大观测值，界限之外的观测值通常单独标为异常点。
par(family = cn_plot_family)
boxplot(
grades,
col = "orange3",
main = "《计算统计》课程成绩箱线图",
ylab = "成绩"
)
散点图用于展示两个定量变量之间的关系。
par(family = cn_plot_family)
plot(
women$height, women$weight,
col = "#2C7FB8",
pch = 16,
xlab = "身高",
ylab = "体重",
main = "身高与体重"
)
ggplot2 与图形语法
ggplot2 基于图形语法，把一幅图看成数据、变量映射、几何对象、坐标系统、标度和主题等组成部分的叠加；关于 ggplot2 的系统介绍可参见 Wickham (2016)。与基础绘图相比，ggplot2 更适合构造层次清楚、便于修改的统计图形。
library(ggplot2)
set.seed(123)
diamonds_small <- diamonds[sample(seq_len(nrow(diamonds)), 1000), ]
ggplot(diamonds_small, aes(x = carat, y = price, color = cut)) +
geom_point(alpha = 0.35) +
geom_smooth(se = FALSE) +
labs(
title = "钻石重量与价格",
x = "重量（carat）",
y = "价格",
color = "切工"
) +
theme_minimal()
这段代码体现了 ggplot2 的基本思路：ggplot() 提供数据和变量映射，geom_point() 加入散点层，geom_smooth() 加入平滑曲线层，labs() 修改标签，theme_minimal() 修改主题。后续章节中，我们将用图形检查算法收敛、模拟误差、残差模式和后验样本。
2.11 进一步阅读
R 语言的历史和设计可参见 Ihaka and Gentleman (1996) 与 Chambers (2008)。R 语言基础可参考 Adler (2010)。可复现研究的基本原则可参见 Gentleman and Temple Lang (2007) 和 Roger D. Peng (2011)；动态文档和 R Markdown 的系统介绍可参见 Xie (2015) 和 Xie, Allaire, and Grolemund (2018)。统计图形和 ggplot2 可参见 Wickham (2016)。
2.12 本章小结
本章介绍了计算统计所需的 R 语言基础。R 中的向量、矩阵、列表和数据框分别适合表示样本、模型矩阵、复杂输出和观测数据；子集提取和向量化运算是 R 编程的核心，但也容易带来隐蔽的维度错误；函数化编程可以把统计计算步骤变成可测试、可复用的模块；随机数种子、相对路径、动态报告和环境记录共同构成可复现计算的基础。后续章节将在这些工具之上讨论更具体的数值方法和统计计算算法。
2.13 思考题
- 为什么说“代码可以运行”并不等于“分析可以复现”？请至少从数据、随机数、软件环境和报告生成四个方面说明。
- 比较向量、矩阵、列表和数据框在统计建模中的作用。模型矩阵、回归系数、模型拟合结果和原始调查数据分别更适合用哪类对象表示？
- 为什么
X[, 1]与X[, 1, drop = FALSE]可能导致不同的后续计算结果？请结合矩阵运算举例说明。 - 在 Monte Carlo 实验中，为什么不应在每一次重复实验中使用同一个随机数种子？
- 函数
normal_loglik()为什么使用 \(\log\sigma\) 作为参数，而不是直接使用 \(\sigma\)？ - 条形图、直方图、箱线图和散点图分别适合回答什么类型的数据问题？
2.14 上机实验（Lab）
- Lab 1：项目结构与动态报告。 建立一个 R 项目目录，包含
data/、R/和output/三个子目录。用 R Markdown 写一份简短报告，读取一个 CSV 文件，计算描述统计量，生成一幅图，并在报告末尾输出sessionInfo()。 - Lab 2：函数化似然计算。 修改本章的
normal_loglik()函数，使其接收参数向量和数据向量，并返回负对数似然。用optim()对模拟正态样本进行极大似然估计，比较估计值与样本均值、样本标准差之间的关系。 - Lab 3：随机数种子与 Monte Carlo 重复。 估计 \(\operatorname{E}(X^2)\)，其中 \(X\sim N(0,1)\)。分别比较以下两种做法：每次重复实验都设置同一个种子；只在实验开始前设置一次种子。解释两种结果的差异。
- Lab 4：图形诊断。 生成一组含有异常值的线性回归模拟数据。分别画散点图、残差图和箱线图，说明每幅图能揭示什么问题。实验报告应包含代码、图形和文字解释。