第 9 章 随机数生成
前面几章主要讨论确定性算法：给定同样的输入，算法每次都会给出同样的输出。从本章开始，我们进入随机模拟方法。随机模拟的基本思想是：如果某个统计量难以直接计算，就通过大量随机样本来近似它。因此，随机数生成是 Monte Carlo、Bootstrap 和 MCMC 的共同起点。
在经济统计中，随机模拟有很多自然用途。例如，我们可以模拟家庭收入分布来评估贫困率估计的不确定性；模拟抽样调查过程来比较不同样本量的标准误；模拟企业违约事件来研究风险模型的误差；也可以在课程项目中用模拟数据检查估计算法是否按理论预期工作。
本章的学习重点有三个。第一，理解计算机中的随机数为什么可以复现；第二，掌握从均匀随机数生成常见分布随机数的基本方法；第三，能够写出清楚的数据生成过程，为后续 Monte Carlo、Bootstrap 和MCMC 章节做准备。
9.1 伪随机数与可重复性
计算机生成的随机数通常不是“真正随机”的，而是由确定性算法产生的伪随机数。给定同一个随机种子，算法会生成完全相同的随机数序列。这一点对统计计算非常重要：模拟结果既要像随机实验，又要能够被教师、同学或审稿人复现。
set.seed(2026)
runif(5)
#> [1] 0.6987 0.5565 0.1401 0.2857 0.5554
set.seed(2026)
runif(5)
#> [1] 0.6987 0.5565 0.1401 0.2857 0.5554两次结果相同，是因为我们把随机数生成器重置到了同一个初始状态。写模拟代码时，通常应在主要模拟开始前明确设置 set.seed()。但也不要在循环内部反复设置同一个种子，否则会不断重复同一批随机数。
R 中可以查看当前随机数生成器类型：
本科阶段不需要深入研究随机数生成器的数学构造，但要形成两个基本习惯：第一，模拟实验要设置随机种子；第二，报告模拟结果时要说明模拟设定，而不是只给出一个无法复现的数字。
一个最简单的伪随机数生成思想是线性同余法：
\[ x_{t+1}=(a x_t+c)\bmod m, \]
其中 \(x_t\) 是第 \(t\) 步的整数状态，\(a\)、\(c\) 和 \(m\) 是算法参数，\(\bmod\) 表示取余。把
\[ u_t=\frac{x_t}{m} \]
看作 \((0,1)\) 上的数，就得到一串伪随机数。下面用很小的参数演示这一点。这个例子只是为了说明机制，实际模拟中不要使用这样简陋的生成器。
lcg <- function(seed, a = 17, c = 43, m = 100, n = 12) {
x <- numeric(n)
x[1] <- seed
for (i in 2:n) {
x[i] <- (a * x[i - 1] + c) %% m
}
x / m
}
lcg(seed = 7)
#> [1] 0.07 0.62 0.97 0.92 0.07 0.62 0.97 0.92 0.07 0.62 0.97 0.92
lcg(seed = 7)
#> [1] 0.07 0.62 0.97 0.92 0.07 0.62 0.97 0.92 0.07 0.62 0.97 0.92两次结果相同，说明“随机种子”本质上就是随机数生成器的初始状态。好的随机数生成器会使用更复杂的状态更新规则，使序列在统计性质上接近真正随机。R 默认的随机数生成器已经足够支持本书中的教学模拟。
还要注意，不应在循环内部反复设置同一个种子。下面的代码会生成 5 个完全相同的数：
bad <- numeric(5)
for (i in 1:5) {
set.seed(99)
bad[i] <- runif(1)
}
bad
#> [1] 0.5847 0.5847 0.5847 0.5847 0.5847正确做法是在模拟开始前设置一次种子，然后让随机数生成器自然向前推进。
9.2 均匀随机数
最基本的随机数是均匀分布随机数。若
\[ U\sim U(0,1), \]
则 \(U\) 在区间 \((0,1)\) 上各处取值机会相同。R 中使用 runif() 生成均匀随机数：
u <- runif(10)
u
#> [1] 0.025131 0.466231 0.861011 0.252501 0.580806 0.005933 0.691866 0.231125
#> [9] 0.848532 0.153835很多其他分布的随机数都可以从均匀随机数变换得到。直观地说，均匀随机数提供了“随机性来源”，而不同的变换把这种随机性塑造成指数分布、正态分布、二项分布等不同形状。
还有一个与逆变换法配套的重要事实：如果 \(X\) 的分布函数为 \(F\)，且 \(F\) 连续，则
\[ F(X)\sim U(0,1). \]
这称为概率积分变换（probability integral transform）。它说明分布函数可以把任意连续随机变量重新标定到 \((0,1)\) 区间上；反过来，分位数函数又可以把均匀随机数变成目标分布随机数。下面用正态随机数做一个简单检查。
set.seed(2027)
x_norm <- rnorm(10000)
u_from_norm <- pnorm(x_norm)
c(mean = mean(u_from_norm),
variance = var(u_from_norm),
theoretical_mean = 1 / 2,
theoretical_variance = 1 / 12)
#> mean variance theoretical_mean
#> 0.49676 0.08335 0.50000
#> theoretical_variance
#> 0.08333如果画出 u_from_norm 的直方图，应接近一条高度为 1 的水平线。这个例子也帮助我们理解 R 中p 函数和 q 函数的关系：pnorm() 把正态变量变成均匀变量，qnorm() 则把均匀变量变回正态变量。
R 对常见分布函数采用统一命名规则。以前缀为例：
| 前缀 | 含义 | 正态分布例子 | Poisson 分布例子 |
|---|---|---|---|
| d | 密度或概率质量函数 | dnorm() | dpois() |
| p | 分布函数 | pnorm() | ppois() |
| q | 分位数函数 | qnorm() | qpois() |
| r | 随机数生成 | rnorm() | rpois() |
例如，rnorm(100) 表示生成 100 个正态随机数，qnorm(0.975) 表示标准正态分布的 97.5% 分位数。熟悉这套命名规则，可以显著减少查函数帮助的成本。
9.3 逆变换法
设随机变量 \(X\) 的分布函数为 \(F(x)\)。如果 \(F\) 连续且严格单调，则
\[ X=F^{-1}(U),\quad U\sim U(0,1) \]
服从分布函数 \(F\)。这是因为
\[ P\{F^{-1}(U)\leq x\}=P\{U\leq F(x)\}=F(x). \]
其中 \(F^{-1}\) 是分位数函数。这个方法称为逆变换法。以指数分布为例，若\(X\sim \operatorname{Exp}(\lambda)\)，则
\[ F(x)=1-\exp(-\lambda x),\quad x\geq 0. \]
令 \(U=F(X)\)，可得
\[ X=-\frac{1}{\lambda}\log(1-U). \]
set.seed(1)
lambda <- 2
u <- runif(10000)
x <- -log(1 - u) / lambda
c(sample_mean = mean(x), theoretical_mean = 1 / lambda)
#> sample_mean theoretical_mean
#> 0.503 0.500R 已经提供了 rexp()，但手写一次逆变换有助于理解随机数生成背后的原理。
逆变换法也可以用于离散分布。若随机变量 \(X\) 取值为 \(x_1,\ldots,x_K\)，对应概率为\(p_1,\ldots,p_K\)，则可以先计算累积概率
\[ P_k=\sum_{j=1}^k p_j, \]
再生成 \(U\sim U(0,1)\)，取满足 \(U\leq P_k\) 的第一个类别。下面手写一个离散抽样函数。
sample_discrete_inverse <- function(values, prob, n) {
prob <- prob / sum(prob)
cum_prob <- cumsum(prob)
u <- runif(n)
values[findInterval(u, c(0, cum_prob), rightmost.closed = TRUE)]
}
set.seed(11)
labor_status <- sample_discrete_inverse(
c("就业", "失业", "退出劳动力市场"),
prob = c(0.72, 0.08, 0.20),
n = 10
)
labor_tab <- data.frame(index = seq_along(labor_status),
status = labor_status)
names(labor_tab) <- c("序号", "劳动力状态")
knitr::kable(labor_tab, row.names = FALSE)| 序号 | 劳动力状态 |
|---|---|
| 1 | 就业 |
| 2 | 就业 |
| 3 | 就业 |
| 4 | 就业 |
| 5 | 就业 |
| 6 | 退出劳动力市场 |
| 7 | 就业 |
| 8 | 就业 |
| 9 | 退出劳动力市场 |
| 10 | 就业 |
R 的 sample() 已经封装了这类操作，但手写版本有助于理解类别抽样如何从均匀随机数出发。
9.4 变换法与正态随机数
正态分布随机数的生成比指数分布更复杂。一个经典方法是 Box–Muller 变换。若\(U_1,U_2\) 独立同分布于 \(U(0,1)\)，则
\[ Z_1=\sqrt{-2\log U_1}\cos(2\pi U_2) \]
服从标准正态分布。这里 \(U_1\) 控制半径，\(U_2\) 控制角度，变换后得到正态形状。
set.seed(2)
n <- 10000
u1 <- runif(n)
u2 <- runif(n)
z <- sqrt(-2 * log(u1)) * cos(2 * pi * u2)
c(mean = mean(z), sd = sd(z))
#> mean sd
#> 0.02117 0.99594实际使用中，我们通常直接调用 rnorm()。学习 Box–Muller 变换的目的不是替代 R 的实现，而是理解：复杂分布的随机数往往由简单随机数经过数学变换得到。
9.5 拒绝抽样简介
有些分布既没有简单的分位数函数，也不容易通过简单变换生成随机数。这时可以使用拒绝抽样。设目标密度为 \(f(x)\)，选择一个容易抽样的建议密度 \(g(x)\)，并找到常数 \(M\)，使得
\[ f(x)\leq M g(x) \]
对所有 \(x\) 成立。拒绝抽样的步骤如下。
拒绝抽样
- 从建议分布 \(g(x)\) 中生成候选值 \(X^\ast\)。
- 生成 \(U\sim U(0,1)\)。
- 若
\[ U\leq \frac{f(X^\ast)}{M g(X^\ast)}, \]
则接受 \(X^\ast\)；否则拒绝并重新抽样。
直观地说，\(M g(x)\) 是罩住目标密度 \(f(x)\) 的“外壳”。在外壳下面按比例接受候选值，被接受的样本就服从目标分布。一个合适的建议分布至少要满足两点：
- 目标分布可能取值的地方，建议分布也必须有正密度。否则某些区域永远抽不到。
- 比值 \(f(x)/g(x)\) 必须有有限上界，也就是存在有限的 \(M\)。实践中常选择尾部比目标分布更厚的建议分布，以避免尾部位置的比值爆炸。
下面用自由度为 3 的 \(t\) 分布作为建议分布，生成标准正态随机数。这个例子主要用于说明算法，实际工作中当然直接使用 rnorm()。
target <- function(x) dnorm(x)
proposal_density <- function(x) dt(x, df = 3)
ratio <- function(x) target(x) / proposal_density(x)
M <- optimize(ratio, interval = c(-10, 10), maximum = TRUE)$objective
M
#> [1] 1.17
1 / M
#> [1] 0.8544这里 \(1/M\) 是理论接受率。由于 \(t_3\) 分布比标准正态分布尾部更厚，它可以作为覆盖正态密度的建议分布；但它也会生成较多极端候选值，这些候选值往往更容易被拒绝。
set.seed(12)
proposal_sample <- function(n) rt(n, df = 3)
n_accept <- 2000
x_accept <- numeric(n_accept)
k <- 0
trial <- 0
while (k < n_accept) {
trial <- trial + 1
x_star <- proposal_sample(1)
alpha <- target(x_star) / (M * proposal_density(x_star))
if (runif(1) <= alpha) {
k <- k + 1
x_accept[k] <- x_star
}
}
c(sample_mean = mean(x_accept),
sample_sd = sd(x_accept),
acceptance_rate = n_accept / trial)
#> sample_mean sample_sd acceptance_rate
#> 0.01658 1.00445 0.84998为什么被接受的样本服从目标分布？直观解释是：候选值先按 \(g(x)\) 的形状出现，再按\(f(x)/(M g(x))\) 的比例保留。因此某个小区间附近被保留下来的概率与
\[ g(x)\frac{f(x)}{M g(x)}=\frac{f(x)}{M} \]
成正比，归一化以后就是目标密度 \(f(x)\)。同时，单个候选值被接受的总概率为
\[ \int \frac{f(x)}{M g(x)}g(x)\,dx =\frac{1}{M}. \]
拒绝抽样的效率取决于建议分布是否接近目标分布。若 \(M\) 很大，接受率就低；若目标分布维度很高，找到合适的建议分布也会变得困难。这正是后续 MCMC 方法的重要动机之一。
从 Monte Carlo 计算角度看，拒绝抽样也留下了一个自然问题：被拒绝的候选值虽然比例不合适，但并非完全没有信息。下一章的重要性抽样会沿着这个思路继续前进：不再简单丢弃候选值，而是用权重修正它们对目标期望的贡献。
9.6 离散分布与抽样
经济统计中经常需要从离散类别中抽样。例如，模拟家庭消费主项类别、企业行业类别或劳动力状态。R 中的 sample() 可以按指定概率抽样。
set.seed(3)
categories <- c("食品", "居住", "交通", "教育", "医疗")
prob <- c(0.30, 0.28, 0.16, 0.14, 0.12)
draw20 <- sample(categories, size = 20, replace = TRUE, prob = prob)
draw20_tab <- data.frame(index = seq_along(draw20),
category = draw20)
names(draw20_tab) <- c("序号", "消费主项")
knitr::kable(
draw20_tab,
row.names = FALSE,
caption = "前 20 个模拟家庭的消费主项"
)| 序号 | 消费主项 |
|---|---|
| 1 | 食品 |
| 2 | 教育 |
| 3 | 居住 |
| 4 | 居住 |
| 5 | 交通 |
| 6 | 交通 |
| 7 | 食品 |
| 8 | 食品 |
| 9 | 居住 |
| 10 | 交通 |
| 11 | 居住 |
| 12 | 居住 |
| 13 | 居住 |
| 14 | 居住 |
| 15 | 教育 |
| 16 | 教育 |
| 17 | 食品 |
| 18 | 交通 |
| 19 | 医疗 |
| 20 | 食品 |
draw1000 <- sample(categories, size = 1000, replace = TRUE, prob = prob)
freq <- prop.table(table(draw1000))
freq_tab <- data.frame(category = names(freq),
proportion = as.numeric(freq))
names(freq_tab) <- c("消费主项", "模拟比例")
knitr::kable(
freq_tab,
row.names = FALSE,
digits = 3,
caption = "1000 次模拟抽样得到的类别比例"
)| 消费主项 | 模拟比例 |
|---|---|
| 交通 | 0.159 |
| 医疗 | 0.117 |
| 居住 | 0.278 |
| 教育 | 0.155 |
| 食品 | 0.291 |
这里 replace = TRUE 表示每次抽样后放回，因此每个家庭的消费主项类别是独立生成的。若设置replace = FALSE，则表示不放回抽样，更接近从一个有限总体中抽取样本。
经济统计中还经常需要模拟有限总体抽样。假设一个城市有 \(N\) 个家庭，我们先生成一个有限总体，再从中不放回抽取样本。样本均值会围绕总体均值波动。
set.seed(13)
N <- 10000
population_income <- rlnorm(N, meanlog = log(8), sdlog = 0.55)
pop_mean <- mean(population_income)
n_sample <- 300
sample_id <- sample(seq_len(N), size = n_sample, replace = FALSE)
sample_mean <- mean(population_income[sample_id])
c(population_mean = pop_mean,
sample_mean = sample_mean,
sampling_error = sample_mean - pop_mean)
#> population_mean sample_mean sampling_error
#> 9.278524 9.271229 -0.007295这个例子区分了两个层次的随机性：有限总体本身是我们模拟出来的；从有限总体中抽取样本又引入了抽样误差。真实抽样调查中，总体不是程序生成的，但抽样误差同样存在。
9.7 一个经济统计模拟数据生成器
为了后续章节使用，我们写一个稍完整的家庭调查模拟函数。它生成地区、收入、就业状态、消费支出和抽样权重。所有数据都是模拟的，不对应任何真实调查。
simulate_households <- function(n) {
region <- sample(c("东部", "中部", "西部"),
size = n, replace = TRUE,
prob = c(0.42, 0.33, 0.25))
region_effect <- ifelse(region == "东部", 0.18,
ifelse(region == "西部", -0.16, 0))
income <- rlnorm(n, meanlog = log(8) + region_effect, sdlog = 0.55)
employed <- rbinom(n, size = 1, prob = plogis(-0.4 + 0.18 * income))
consumption <- 1.2 + 0.62 * income + 0.8 * employed +
ifelse(region == "东部", 0.4, ifelse(region == "西部", -0.25, 0)) +
rnorm(n, sd = 1.5)
consumption <- pmax(consumption, 0.2)
weight <- runif(n, min = 0.7, max = 1.5)
data.frame(
region = region,
income = income,
employed = employed,
consumption = consumption,
weight = weight
)
}
set.seed(4)
dat <- simulate_households(5)
knitr::kable(dat, row.names = FALSE, digits = 4)| region | income | employed | consumption | weight |
|---|---|---|---|---|
| 中部 | 5.620 | 1 | 6.334 | 1.1537 |
| 东部 | 19.763 | 0 | 13.877 | 0.8912 |
| 东部 | 4.308 | 1 | 5.645 | 1.4024 |
| 东部 | 7.019 | 0 | 5.884 | 1.2236 |
| 西部 | 17.228 | 1 | 12.483 | 1.0859 |
这个函数包含了一个数据生成过程。收入服从右偏分布，不同地区有不同收入水平；就业概率随收入上升而提高；消费支出由收入、就业状态、地区和随机扰动共同决定。这样的模拟设定虽然简化，却非常适合用来学习估计、预测和不确定性评估。
模拟数据生成器应尽量写成函数，而不是把代码散落在脚本各处。这样做有三个好处：第一，便于改变样本量；第二，便于在 Monte Carlo 实验中重复调用；第三，便于在报告中清楚说明数据生成过程。
9.8 随机数生成的检查
生成随机数后，应做一些简单检查。检查不是为了证明随机数“绝对正确”，而是为了发现明显的程序错误。常见检查包括样本均值、样本方差、取值范围和类别比例。
set.seed(5)
sim_dat <- simulate_households(1000)
c(mean_income = mean(sim_dat$income),
sd_income = sd(sim_dat$income),
employment_rate = mean(sim_dat$employed),
mean_consumption = mean(sim_dat$consumption))
#> mean_income sd_income employment_rate mean_consumption
#> 9.967 6.199 0.743 8.077
region_prop <- prop.table(table(sim_dat$region))
region_tab <- data.frame(region = names(region_prop),
proportion = as.numeric(region_prop))
names(region_tab) <- c("地区", "比例")
knitr::kable(
region_tab,
row.names = FALSE,
digits = 3
)| 地区 | 比例 |
|---|---|
| 东部 | 0.421 |
| 中部 | 0.320 |
| 西部 | 0.259 |
如果某个变量的均值、比例或范围明显不符合设定，就应回到数据生成函数检查概率、参数和变量变换。模拟研究的可信度首先来自透明、可检查的数据生成过程。