19 随机模拟

19.1 随机数

随机模拟是统计研究的重要方法, 另外许多现代统计计算方法(如MCMC)也是基于随机模拟的。 R中提供了多种不同概率分布的随机数函数, 可以批量地产生随机数。 一些R扩展包利用了随机模拟方法,如boot包进行bootstrap估计。

所谓随机数,实际是“伪随机数”, 是从一组起始值(称为种子), 按照某种递推算法向前递推得到的。 所以,从同一种子出发,得到的随机数序列是相同的。

为了得到可重现的结果, 随机模拟应该从固定不变的种子开始模拟。 用set.seed(k)指定一个编号为k的种子, 这样每次从编号k种子运行相同的模拟程序可以得到相同的结果。

还可以用set.seed()加选项kind=指定后续程序要使用的随机数发生器名称, 用normal.kind=指定要使用的正态分布随机数发生器名称。

R提供了多种分布的随机数函数,如runif(n)产生n个标准均匀分布随机数, rnorm(n)产生n个标准正态分布随机数。 例如:

注意因为没有指定种子,每次运行会得到不同的结果。

在R命令行运行

可以查看R中提供的不同概率分布。

19.2 sample()函数

sample()函数从一个有限集合中无放回或有放回地随机抽取, 产生随机结果。

例如,为了设随机变量\(X\)取值于\(\{\text{正面}, \text{反面} \}\), 且\(P(X=\text{正面}) = 0.7 = 1 - P(X = \text{反面})\), 如下程序产生\(X\)的10个随机抽样值:

sample()的选项size指定抽样个数, prob指定每个值的概率, replace=TRUE说明是有放回抽样。

如果要做无放回等概率的随机抽样, 可以不指定probreplace(缺省是FALSE)。 比如,下面的程序从1:10随机抽取4个:

如果要从\(1:n\)中等概率无放回随机抽样直到每一个都被抽过,只要用如:

这实际上返回了\(1:10\)的一个重排。

dplyr包的sample_n()函数与sample()类似, 但输入为数据框, 输出为随机抽取的数据框行子集。

19.3 随机模拟示例

19.3.1 估计希望值

设随机变量或随机向量\(X\)有复杂的分布, 使得希望\(\theta = E h(X)\)很难计算。 如果可以得到\(X\)\(N\)个独立同分布抽样值\(X_1, X_2, \dots, X_N\), 则可估计\(\theta\)为: \[ \hat\theta = \frac{1}{N} \sum_{i=1}^N h(X_i) \] 这时\(E \hat\theta = \theta\), 估计无偏; 均方误差为 \[\begin{aligned} \text{MSE} =& E|\hat\theta - \theta|^2 = \text{Var}(\hat\theta) = \text{SE}^2 \\ =& \frac{\text{Var}(X)}{N} \approx \frac{S_N^2}{N} \end{aligned}\] 其中 \[ S_N^2 = \frac{1}{N-1} \sum_{i=1}^N (h(X_i) - \hat\theta)^2 \] 是样本方差。 由强大数律可知\(N\to\infty\)\(\hat\theta\)依概率1收敛到\(\theta\), 由中心极限定理可知\(N\)充分大时\(\hat\theta\)有近似正态分布 \[ \text{N}(\theta, \text{SE}^2) . \]

例如,考虑正方形区域\(\Omega = \{(x,y): x \in [0,1], y \in [0,1] \}\), 以及其中的四分之一扇形\(A = \{(x,y): (x,y) \in \Omega, x^2 + y^2 \leq 1 \}\)。 设\(\boldsymbol X\)服从\(\Omega\)上的均匀分布, 令 \[ Y = \begin{cases} 1, & \text{当} \boldsymbol X \in A, \\ 0, & \text{其它} \end{cases} \]\(Y\)服从两点分布, 概率为 \[ p = E Y = \frac{\frac{1}{4} \pi 1^2}{1^2} = \frac{\pi}{4} \]

\(\Omega\)中投入\(N=10000\)个点, 得到\(N\)\(Y_i, i=1,2,\dots,N\)的值,估计\(\pi\)\[ \hat\pi = 4 \bar Y = \frac{4}{N} \sum_{i=1}^N Y_i \] \(E \hat\pi = \pi\), 估计的根均方误差可估计为 \[ \text{RMSE} = \sqrt{E | \hat\pi - \pi |^2} = 4 \frac{\sqrt{\text{Var}(Y)}}{\sqrt{N}} \approx 4 \frac{S_N}{\sqrt{N}} . \]

程序如下:

估计的标准误差(SE)还是比较大, 提高到\(N=10^6\),精度可以增加一位小数:

19.3.2 线性回归模拟

考虑如下线性回归模型 \[\begin{aligned} y = 10 + 2 x + \varepsilon, \ \varepsilon \sim \text{N}(0, 0.5^2) . \end{aligned}\] 假设有样本量\(n=10\)的一组样本, R函数lm()可以 可以得到截距\(a\), 斜率\(b\)的估计\(\hat a, \hat b\), 以及相应的标准误差\(\text{SE}(\hat a)\), \(\text{SE}(\hat b)\)。 样本可以模拟产生。

模型中的自变量\(x\)可以用随机数产生, 比如,用sample()函数从\(1:10\)中随机有放回地抽取\(n\)个。 模型中的随机误差项\(\varepsilon\)可以用rnorm()产生。 产生一组样本的程序如:

如下程序计算线性回归:

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      10.080        1.995

如下程序计算线性回归的多种统计量:

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76774 -0.14896  0.02445  0.33051  0.49779 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.08026    0.22313   45.18 6.37e-11 ***
## x            1.99518    0.04545   43.90 8.00e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4286 on 8 degrees of freedom
## Multiple R-squared:  0.9959, Adjusted R-squared:  0.9953 
## F-statistic:  1927 on 1 and 8 DF,  p-value: 8.004e-11

如下程序返回一个矩阵, 包括\(a, b\)的估计值、标准误差、t检验统计量、检验p值:

##              Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 10.080258 0.22313107 45.17640 6.365149e-11
## x            1.995183 0.04545185 43.89664 8.003683e-11

如下程序把上述矩阵的前两列拉直成一个向量返回:

## [1] 10.08025775  1.99518338  0.22313107  0.04545185

这样得到 \(\hat a, \hat b, \text{SE}(\hat a), \text{SE}(\hat b)\)这四个值。

replicate(, 复合语句)实行多次模拟, 返回向量或矩阵结果, 返回矩阵时,每列是一次模拟的结果。 下面是线性回归整个模拟程序,写成了一个函数。

运行测试:

## 1000 次模拟的平均值:
##         a         b      SE.a      SE.b 
## 9.9970476 1.9994490 0.3639505 0.0592510 
## 1000 次模拟的标准差:
##          a          b       SE.a       SE.b 
## 0.37974881 0.06297733 0.11992515 0.01795926

可以看出,标准误差作为\(\hat a, \hat b\)的标准差估计, 与多次模拟得到多个\(\hat a, \hat b\)样本计算得到的标准差估计是比较接近的。 结果中\(\text{SE}(\hat a)\)的平均值为0.363, 1000次模拟的\(\hat a\)的样本标准差为0.393,比较接近; \(\text{SE}(\hat b)\)的平均值为0.0594, 1000次模拟的\(\hat b\)的样本标准差为0.0637,比较接近。

19.3.3 核密度的bootstrap置信区间

R自带的数据框faithful内保存了美国黄石国家公园Faithful火山的272次爆发持续时间和间歇时间。 为估计爆发持续时间的密度,可以用核密度估计方法, R函数density可以实行此估计, 返回\(N\)个格子点上的密度曲线坐标:

这个密度估计明显呈现出双峰形态。

核密度估计是统计估计, 为了得到其置信区间(给定每个\(x\)坐标,真实密度\(f(x)\)的单点的置信区间), 采用如下非参数bootstrap方法:

重复\(B=10000\)次, 每次从原始样本中有重复地抽取与原来大小相同的一组样本, 对这组样本计算核密度估计, 结果为\((x_i, y_i^{(j)}), i=1,2,\dots,N, j=1,2,\dots, B\), 每组样本估计\(N\)个格子点的密度曲线坐标, 横坐标不随样本改变。

对每个横坐标\(x_i\), 取bootstrap得到的\(B\)\(y_i^{(j)}, j=1,2,\dots,B\)的0.025和0.975样本分位数, 作为真实密度\(f(x_i)\)的bootstrap置信区间。

在R中利用replicate()函数实现:

程序中用set.seed(1)保证每次运行得到的结果是不变的, replicate()函数第一参数是重复模拟次数, 第二参数是复合语句, 这些语句是每次模拟要实行的计算。 在每次模拟中, 用带有replace=TRUE选项的sample()函数从样本中有放回地抽样得到一组bootstrap样本, 每次模拟的结果是在指定格子点上计算的核密度估计的纵坐标。 replicate()的结果为一个矩阵, 每一列是一次模拟得到的纵坐标集合。 对每个横坐标格子点,用quantile()函数计算\(B\)个bootstrap样本的2.5%和97.5%分位数, 循环用apply()函数表示。 polygon()函数指定一个多边形的顺序的顶点坐标用col=指定的颜色填充, 本程序中实现了置信下限与置信上限两条曲线之间的颜色填充。 lines()函数绘制了与原始样本对应的核密度估计曲线。