第 9 章 随机数生成
前面几章主要讨论确定性算法:给定同样的输入,算法每次都会给出同样的输出。从本章开始,我们进入随机模拟方法。随机模拟的基本思想是:如果某个统计量难以直接计算,就通过大量随机样本来近似它。因此,随机数生成是 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 中可以查看当前随机数生成器类型:
本科阶段不需要深入研究随机数生成器的数学构造,但要形成两个基本习惯:第一,模拟实验要设置随机种子;第二,报告模拟结果时要说明模拟设定,而不是只给出一个无法复现的数字。
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很多其他分布的随机数都可以从均匀随机数变换得到。直观地说,均匀随机数提供了“随机性来源”,而不同的变换把这种随机性塑造成指数分布、正态分布、二项分布等不同形状。
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(),但手写一次逆变换有助于理解随机数生成背后的原理。
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 离散分布与抽样
经济统计中经常需要从离散类别中抽样。例如,模拟家庭消费主项类别、企业行业类别或劳动力状态。R 中的 sample() 可以按指定概率抽样。
set.seed(3)
categories <- c("食品", "居住", "交通", "教育", "医疗")
prob <- c(0.30, 0.28, 0.16, 0.14, 0.12)
sample(categories, size = 20, replace = TRUE, prob = prob)
#> [1] "<U+98DF><U+54C1>" "<U+6559><U+80B2>" "<U+5C45><U+4F4F>" "<U+5C45><U+4F4F>"
#> [5] "<U+4EA4><U+901A>" "<U+4EA4><U+901A>" "<U+98DF><U+54C1>" "<U+98DF><U+54C1>"
#> [9] "<U+5C45><U+4F4F>" "<U+4EA4><U+901A>" "<U+5C45><U+4F4F>" "<U+5C45><U+4F4F>"
#> [13] "<U+5C45><U+4F4F>" "<U+5C45><U+4F4F>" "<U+6559><U+80B2>" "<U+6559><U+80B2>"
#> [17] "<U+98DF><U+54C1>" "<U+4EA4><U+901A>" "<U+533B><U+7597>" "<U+98DF><U+54C1>"
table(sample(categories, size = 1000, replace = TRUE, prob = prob)) / 1000
#>
#> <U+4EA4><U+901A> <U+533B><U+7597> <U+5C45><U+4F4F> <U+6559><U+80B2>
#> 0.159 0.117 0.278 0.155
#> <U+98DF><U+54C1>
#> 0.291这里 replace = TRUE 表示每次抽样后放回,因此每个家庭的消费主项类别是独立生成的。若设置replace = FALSE,则表示不放回抽样,更接近从一个有限总体中抽取样本。
9.6 一个经济统计模拟数据生成器
为了后续章节使用,我们写一个简单的家庭调查模拟函数。它生成收入、就业状态和消费支出。所有数据都是模拟的,不对应任何真实调查。
simulate_households <- function(n) {
income <- rlnorm(n, meanlog = log(8), sdlog = 0.55)
employed <- rbinom(n, size = 1, prob = plogis(-0.4 + 0.18 * income))
consumption <- 1.2 + 0.62 * income + 0.8 * employed + rnorm(n, sd = 1.5)
consumption <- pmax(consumption, 0.2)
data.frame(
income = income,
employed = employed,
consumption = consumption
)
}
set.seed(4)
dat <- simulate_households(5)
dat
#> income employed consumption
#> 1 9.013 1 7.419
#> 2 5.936 1 5.999
#> 3 13.060 1 11.165
#> 4 11.103 0 12.146
#> 5 19.669 1 14.157这个函数包含了一个数据生成过程。收入服从右偏分布,就业概率随收入上升而提高,消费支出由收入、就业状态和随机扰动共同决定。这样的模拟设定虽然简化,却非常适合用来学习估计、预测和不确定性评估。