第 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 中可以查看当前随机数生成器类型:

RNGkind()
#> [1] "Mersenne-Twister" "Inversion"        "Rejection"

本科阶段不需要深入研究随机数生成器的数学构造,但要形成两个基本习惯:第一,模拟实验要设置随机种子;第二,报告模拟结果时要说明模拟设定,而不是只给出一个无法复现的数字。

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.500

R 已经提供了 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

这个函数包含了一个数据生成过程。收入服从右偏分布,就业概率随收入上升而提高,消费支出由收入、就业状态和随机扰动共同决定。这样的模拟设定虽然简化,却非常适合用来学习估计、预测和不确定性评估。

9.7 本章小结

随机数生成是随机模拟方法的基础。伪随机数由确定性算法产生,因此可以通过随机种子复现。均匀随机数是许多分布随机数的起点,逆变换法和变换法说明了如何从简单随机性构造复杂分布。对于经济统计教学,模拟数据不仅可以帮助我们理解算法,也能在不依赖外部数据来源的情况下检验统计方法。

9.8 练习

  1. 用逆变换法生成 10000 个参数为 \(\lambda=0.5\) 的指数分布随机数,比较样本均值和理论均值。
  2. 修改 simulate_households() 中收入分布的 sdlog,观察收入右偏程度如何变化。
  3. sample() 模拟 1000 名劳动者的就业状态,设就业、失业、退出劳动力市场的概率分别为 0.7、0.1、0.2。
  4. 思考为什么在模拟报告中只写“随机生成数据”是不够的。