10 随机模拟介绍

在用数学模型， 包括概率统计模型处理实际应用中的问题时， 我们希望建立的模型能够尽可能地符合实际情况。 但是，实际情况是错综复杂的， 如果一味地要求模型与实际完全相符， 会导致模型过于复杂， 以至于不能进行严格理论分析， 结果导致模型不能使用。 所以，实际建模时会忽略许多细节， 增加一些可能很难验证的理论假设， 使得模型比较简单，可以用数学理论进行分析研究。

这样，简化的模型就可以与实际情况有较大的差距， 即使我们对模型进行了完美的理论分析， 也不能保证分析结果是可信的。 这一困难可以用随机模拟的方法解决。

模拟是指把某一现实的或抽象的系统的某种特征或部分状态， 用另一系统（称为模拟模型）来代替或近似。 为了解决某问题， 把它变成一个概率模型的求解问题， 然后产生符合模型的大量随机数， 对产生的随机数进行分析从而求解问题， 这种方法叫做随机模拟方法， 又称为蒙特卡洛(Monte Carlo)方法。

例如，一个交通路口需要找到一种最优的控制红绿灯信号的办法， 使得通过路口的汽车耽搁的平均时间最短， 而行人等候过路的时间不超过某一给定的心理极限值。 十字路口的信号共有四个方向， 每个方向又分直行、左转、右转。 因为汽车和行人的到来是随机的， 我们要用随机过程来描述四个方向的汽车到来和路口的行人到来过程。 理论建模分析很难解决这个最优化问题。 但是， 我们可以采集汽车和行人到来的频率， 用随机模拟方法模拟汽车和行人到来的过程， 并模拟各种控制方案，记录不同方案造成的等待时间， 通过统计比较找出最优的控制方案。

随机模拟中的随机性可能来自模型本身的随机变量， 比如上面描述的汽车和行人到来， 也可能是把非随机的问题转换为概率模型的特征量估计问题从而用随机模拟方法解决。

例10.1 (用随机模拟估计圆周率) 为了计算圆周率\(\pi\)的近似值可以使用随机模拟方法。

如果向正方形\(D = \{(x,y): x \in [-1,1], y \in [-1,1]\)内随机等可能投点， 落入单位圆\(C = \{(x,y): x^2 + y^2 \leq 1 \}\)的概率为面积之比 \(p = \frac{\pi}{4}\)。 如果独立重复地投了\(N\)个点， 落入\(C\)中的点的个数\(\xi\)的平均值为\(E \xi = p N\)， 由概率的频率解释， \[\begin{aligned} \frac{\xi}{N} \approx& \frac{\pi}{4}, \quad \pi \approx \hat \pi = \frac{4 \xi}{N} \end{aligned}\] 可以这样给出\(\pi\)的近似值。

R程序实现:

sim.pi <- function(N=10000){
  x <- runif(N, -1, 1)
  y <- runif(N, -1, 1)
  p <- mean((x^2 + y^2) <= 1)
  pi.hat <- 4*p
  pi.hat
}

测试:

set.seed(101)
sim.pi(1000000)
## [1] 3.145352
sim.pi(1000000)
## [1] 3.139572

上述模拟的Julia语言版本：

using Random, Distributions
function sim_pi(N=10000)
    x = rand(Uniform(-1.0,1.0), N)
    y = rand(Uniform(-1.0,1.0), N)
    p = mean((x .^ 2 + y .^ 2) .<= 1)
    pihat = 4*p
    pihat
end

## 测试：
Random.seed!(101)
sim_pi(1_000_000)
## 3.140552
sim_pi(1_000_000)
## 3.14062

※※※※※

随机模拟方法会引入所谓随机模拟误差： 上例中估计的\(\hat\pi\)实际是随机的， 如果再次独立重复投\(N\)个点， 得到的\(\hat\pi\)和上一次结果会有不同。 这是随机模拟方法的特点， 即结果具有随机性。 因为结果的随机性导致的误差叫做随机模拟误差。

使用随机模拟方法，我们必须了解随机模拟误差的大小， 这样我们才能够设计合适的重复试验次数来控制随机模拟误差。 比如，这个例子中\(\xi\)服从B(\(N, \frac{\pi}{4}\))分布， 有 \[\begin{aligned} \text{Var}(\hat\pi) =& \frac{\pi(4-\pi)}{N}, \end{aligned}\] 由中心极限定理，\(\hat\pi\)近似服从N(\(\pi\), \(\frac{\pi(4-\pi)}{N}\))分布， 所以随机模拟误差的幅度大约在\(\pm 2 \sqrt{\frac{\pi(4-\pi)}{N}}\) (随机模拟误差95%以上落入此区间)。

例10.2 用上面的例子计算随机模拟误差，并计算95%的误差界限绝对值。

NN <- 10^(2:6)
xmat <- matrix(0, length(NN), 3)
xmat[,1] <- NN
colnames(xmat) <- c("N", "随机模拟误差", "95%误差界")
for(i in seq(nrow(xmat))){
  xmat[i,"随机模拟误差"] <- sim.pi(NN[i]) - pi
  xmat[i, "95%误差界"] <- 2*sqrt(pi*(4-pi)/NN[i])
}
xmat
##          N 随机模拟误差   95%误差界
## [1,] 1e+02  0.258407346 0.328436674
## [2,] 1e+03 -0.013592654 0.103860796
## [3,] 1e+04  0.009207346 0.032843667
## [4,] 1e+05 -0.005552654 0.010386080
## [5,] 1e+06  0.001411346 0.003284367

Julia语言版本：

Random.seed!(101)
NN = 10 .^ (2:6)
xmat = zeros(length(NN), 3)
xmat[:,1] = NN
for i=eachindex(NN)
  xmat[i,2] = sim_pi(NN[i]) - pi
  xmat[i,3] = 2*sqrt(pi*(4-pi)/NN[i])
end
xmat
## 结果略去。

※※※※※

这样的随机误差幅度也是随机模拟误差的典型情况， 样本量增大100倍，误差界才减小到原来的\(\frac{1}{10}\)。

随机模拟方法虽然避免了复杂的理论分析， 但是其结果具有随机性，精度很难提高： 为了增加一位小数点的精度，即误差减小到原来的\(\frac{1}{10}\)， 重复试验次数需要增加到原来的100倍。 随机模拟方法有如下特点:

  • 应用面广，适应性强。只要问题能够清楚地描述出来， 就可以编写模拟程序产生大量数据， 通过分析模拟数据解决问题。

  • 算法简单, 容易改变条件, 但计算量大 。

  • 结果具有随机性，精度较低(一般为\(O(\frac{1}{\sqrt{N}})\)级)。

  • 模拟结果的收敛服从概率论规律。

  • 对维数增加不敏感。在计算定积分时， 如果使用传统数值算法， 维数增加会造成计算时间指数增加， 但是如果使用随机模拟方法计算， 则维数增加常常仅仅使计算时间随维数\(d\)增加而线性地增加。

随机模拟用在科学研究中， 常常作为探索性试验来使用。 假设科学家有了一个新的模型或技术的想法， 但是不知道它的效果怎样所以还没有对其进行深入的理论分析， 就可以用随机模拟方法大量地重复生成模拟数据， 根据多次重复的总体效果来判断这种模型或技术的性能。 如果模拟获得了好的结果， 再进行深入理论分析对模型进行完善； 如果模拟发现了这个模型的缺点， 可以进行有针对性的修改， 或者考虑转而其它解决办法。

随机模拟在科学研究中的另一种作用是说明新的模型或技术的有效性。 在公开发表的统计学论文中，已经有一半以上的文章包括随机模拟结果 （也叫数值结果）， 用来辅助说明自己提出的模型或方法的有效性。 有时因为对模型或方法很难进行彻底的理论分析， 仅仅使用大量的随机模拟结果来说明模型或方法的有效性。 当然，因为模型都是有可变参数的， 随机模拟只能针对某些参数组合给出结果， 所以，一般认为仅有模拟结果而没有理论分析结果的研究论文是不全面的。 第16章给出一个用随机模拟说明统计技术优良性的例子。

除了以上应用， 随机模拟还是许多新的统计方法的主要工具， 例如， bootstrap置信区间和bootstrap偏差修正，置换检验，MCMC。 利用大量计算机计算（包括随机模拟）来进行统计推断的统计学分支叫做 “计算统计”（computational statistics）， 在本章后面各节将介绍随机模拟的一些应用和技巧。