40 随机模拟

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

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

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

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

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

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

在R命令行运行

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

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

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

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

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

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

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

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

设随机变量或随机向量\(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}(h(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}} . \]

程序如下：