15 Simulation, bootstrap statistical methods, and permutation tests
15.1 引言
在本章中，我们将介绍两种强大的现代统计技术：自助法（bootstrap）统计方法和置换（permutation）检验。这两种方法都是非参数方法，他们不对任何潜在概率分布做特定假设。自助法使我们能够衡量参数估计量的有效性，而置换检验则提供了检验某些统计假设的新方法。然而，这两种方法在实施过程中都需要大量的计算。执行所需计算的最有效且高效的方法是使用 模拟（simulation），这也是本章的第三个主题。
在 小节 15.2 中，我们将介绍模拟的关键——随机数（random numbers），我们将展示如何使用随机数来生成随机排列和随机子集。在 小节 15.2.1 中，我们将介绍用于计算近似期望的蒙特卡罗（Monte Carlo）模拟方法。在 ?sec-15_3 中，我们将介绍自助法统计，并展示如何应用蒙特卡罗模拟方法来进行所需的分析。在 ?sec-15_4 中，我们将讨论置换检验：一种用于确定数据序列是否来自单一总体分布的非参数检验方法。在后续的各节中，我们将继续研究模拟。在 ?sec-15_5 和 ?sec-15_6 中，我们将展示如何使用随机数来生成任意连续、离散随机变量的值，在 ?sec-15_7 中，我们将研究何时可以结束蒙特卡罗模拟过程。
15.2 随机数
我们称 \((0, 1)\)-均匀分布随机变量的值为随机数（random number）。过去，我们常常使用机械设备来生成随机数，现在，我们通常使用随机数生成器来生成一个伪随机数序列。随机数生成器首先需要一个初始值 \(x_0\)——种子（seed），然后通过指定的正整数 \(a\)、\(c\) 和 \(m\)，并按照以下的递归方式来确定伪随机序列的值：
\[ x_{n + 1}=(ax_n + c)\bmod m, \quad n\geq0 \]
上述表达式意味着 \(x_{n + 1}\) 是 \(ax_n + c\) 除以 \(m\) 所得的余数。因此，每个 \(x_n\) 是 \(0, 1, \ldots, m - 1\) 中的一个值，并且把 \(x_n/m\) 当作随机数。可以证明，对于合适的 \(a\)、\(c\) 和 \(m\)，上述方法会产生一个序列，该序列看起来就好像是通过观察独立的 \((0, 1)\)-均匀分布随机变量的值而生成的一样。出于这个原因，我们称 \(x_n/m\)（\(n\geq1\)）为伪随机数（pseudo random numbers）。
例 15.1 如果 \(a = 3\)，\(c=7\)，\(m=23\)，当 \(x_0 = 2\) 时，
\[ \begin{align} x_1 &= 3(2) + 7 \quad mod \quad 23 = 13 \\ x_2 &= 3(13) + 7 \quad mod \quad 23 = 0 \\ x_3 &= 3(0) + 7 \quad mod \quad 23 = 7 \\ x_4 &= 3(7) + 7 \quad mod \quad 23 = 5 \\ x_5 &= 3(5) + 7 \quad mod \quad 23 = 22 \end{align} \]
因此，当使用 \(x_0 = 2\) 作为种子时，我们得到了如下的伪随机数序列：\(13/23, 0, 7/23, 5/23, 22/23, \cdots\) \(\blacksquare\)
大多数计算机都有内置的随机数生成器，我们可以使用这些内置的随机数生成器生成伪随机数，并且生成伪随机数是 模拟 过程的起点。此外，我们会假定这些伪随机数实际上就是真正的随机数。也就是说，我们会假定这个随机数序列实际上是来自 \((0, 1)\)-均匀分布的一个样本值序列。
对于任何 模拟 研究而言，随机数都是其关键。我们将在接下来的 例 15.2 中展示随机数的重要性。
例 15.2 假设我们想要生成数字 \(1,2,\ldots,n\) 的一个排列，使得所有 \(n!\) 种可能的排列都具有相同的可能性。我们可以：
- 先从数字 \(1,2,\ldots,n\) 中随机选择一个数字，并将该数字放在位置 \(n\)
- 然后从剩下的 \(n - 1\) 个数字中随机选择一个数字，并将其放在位置 \(n - 1\)
- 接着从剩下的 \(n - 2\) 个数字中随机选择一个数字，并将其放在位置 \(n - 2\)
- 以此类推，直到生成整个排列序列（这里的“随机选择”意味着每个可能的选择都有相同的可能性）
然而，我们无需考虑还有哪些确切的数字需要选择，一种方便且有效的方式是：将所有的数字存储在一个有序列表中，然后随机选择有序序列的位置而不是数字本身。也就是说，从数字 \(1,2,\ldots,n\) 的任何排列 \(r_1,r_2,\ldots,r_n\) 开始：
- 随机选择位置 \(1,\ldots,n\) 中的一个，并将该位置的数字与位置 \(n\) 的数字互换
- 然后随机选择位置 \(1,\ldots,n - 1\) 中的一个，并将该位置的数字与位置 \(n - 1\) 的数字互换
- 接着随机选择位置 \(1,\ldots,n - 2\) 中的一个，并将该位置的数字与位置 \(n - 2\) 的数字互换
- 以此类推
为了实现上述操作，我们需要生成一个随机变量，该随机变量能够等可能地取 \(1,\ldots,k\) 中的任何一个值。为了做到这一点，令 \(U\) 表示一个随机数——即 \(U\) 在 \((0,1)\) 上均匀分布——并令 \(\text{Int}(kU)\) 为 \(kU\) 的整数部分——也就是说，\(\text{Int}(kU)\) 是小于或等于 \(kU\) 的最大整数。那么，对于 \(i = 1,\ldots,k\) ：
\[ \begin{align*} P[\text{Int}(kU)+1 = i]&=P[\text{Int}(kU)=i - 1]\\ &=P(i - 1\leq kU\lt i)\\ &=P(\frac{i - 1}{k}\leq U\lt\frac{i}{k})\\ &=1/k \end{align*} \]
因此， \(\text{Int}(kU)+1\) 等可能地取 \(1,\ldots,k\) 中的任何一个值。
因此，生成数字 \(1,2,\ldots,n\) 的随机排列的算法可以写成如下形式：
- 令 \(r_1,r_2,\ldots,r_n\) 为数字 \(1,2,\ldots,n\) 的任何一个排列。（例如，我们可以令 \(r_j = j\)，\(j = 1,\ldots,n\)。）
- 令 \(k = n\)。（要放在位置 \(k\) 的数字有待确定。）
- 生成一个随机数 \(U\)，并令 \(I=\text{Int}(kU)+1\)。
- 互换 \(r_I\) 和 \(r_k\) 的值。
- 令 \(k = k - 1\)。
- 如果 \(k\gt1\)，转到步骤 3；如果 \(k = 1\)，转到步骤 7。
- \(r_1,\ldots,r_n\) 就是所需的排列。
例如，假设 \(n = 4\)，初始排列是 \(1, 2, 3, 4\)。如果 \(I\) 的第一个值——等可能地为 \(1, 2, 3, 4\) 中的任何一个数——是 \(2\)，那么位置 \(2\) 的数字与位置 \(4\) 的数字互换，得到新的排列 \(1, 4, 3, 2\)。如果 \(I\) 的下一个值——等可能地为 \(1, 2, 3\) 中的任何一个数——是 \(3\)，那么位置 \(3\) 的数字与位置 \(3\) 的数字互换（即数字不变），所以排列仍然是 \(1, 4, 3, 2\)。如果 \(I\) 的最后一个值——等可能地为 \(1, 2\) 中的任何一个数——是 \(1\)，那么位置 \(1\) 的数字与位置 \(2\) 的数字互换，得到最终排列 \(4, 1, 3, 2\)。
上述算法的一个重要性质是，可用于从集合 \(\{1, 2, \ldots, n\}\) 中生成大小为 \(r\) 的随机子集。对于 \(r\leq n/2\)，只需按照上述算法进行操作，直到确定了最后 \(r\) 个位置（即位置 \(n, n - 1, \ldots, n - r + 1\)）中的元素，然后将这些位置中的数字作为大小为 \(r\) 的随机子集。对于 \(r > n/2\)，与其直接选择 \(r\) 个数构成子集，选择不在子集中的 \(n - r\) 个数会更快。所以在这种情况下，按照上述算法进行操作，直到确定了最后 \(n - r\) 个位置的数字，然后将剩下的数字作为大小为 \(r\) 的随机子集。\(\blacksquare\)
在 C++ 标准模板库（STL）中，std::shuffle() 是依赖于随机数生成器的用于随机打乱（洗牌）一个范围内元素顺序的算法。std::shuffle() 的具体算法实现基于 Fisher-Yates 洗牌算法（Knuth 洗牌算法），其核心思想是：从数组末尾开始，逐个与前面的随机位置进行交换，确保最终的排列是完全随机的。具体的代码思路如下所示：
template <typename RandomIt, typename URBG>
void shuffle(RandomIt first, RandomIt last, URBG&& g) {
using diff_t = typename std::iterator_traits<RandomIt>::difference_type;
diff_t n = last - first;
for (diff_t i = n - 1; i > 0; --i) {
std::uniform_int_distribution<diff_t> dist(0, i); // 生成 [0, i] 范围内的随机数
diff_t j = dist(g); // 使用提供的随机数引擎 g
std::swap(first[i], first[j]); // 交换元素
}
}15.2.1 蒙特卡罗模拟方法
假设我们想要计算统计量 \(h(X_1, X_2, \ldots, X_n)\) 的期望，其中 \(X_1, X_2, \ldots, X_n\) 是密度函数为 \(f(x)\) 的独立同分布随机变量。根据 \(X_1, X_2, \ldots, X_n\) 的联合密度函数：
\[ f(x_1, \ldots, x_n) = f(x_1)f(x_2)\cdots f(x_n) \]
我们得到：
\[ \begin{align*} E&[h(X_1, X_2, \ldots, X_n)] \\ &=\int\int\cdots\int h(x_1, \ldots, x_n)f(x_1)f(x_2)\cdots f(x_n)dx_1dx_2\cdots dx_n \end{align*} \]
对于上式的计算而言，其困难在于通常无法通过解析计算（analytically compute）的方法来计算多重积分，同时在指定精度内对上式进行数值计算也很困难。但是，我们仍然可以通过模拟（simulation）来近似计算 \(E[h(X_1, X_2, \ldots, X_n)]\)。
为了实现模拟近似，我们首先生成 \(n\) 个独立随机变量 \(X_1^1, X_2^1, \ldots, X_n^1\) 的值，其中每个随机变量的密度函数都是 \(f\)，然后计算
\[ Y_1 = h(X_1^1, X_2^1, \ldots, X_n^1) \]
接着生成第二组随机变量的值，该组中同样包含 \(n\) 个密度函数为 \(f\) 且与第一组中的随机变量独立的随机变量的值。第二组随机变量记作 \(X_1^2, X_2^2, \ldots, X_n^2\)，然后计算
\[ Y_2 = h(X_1^2, X_2^2, \ldots, X_n^2) \]
重复执行上述操作，直到生成了 \(r\) 组密度函数为 \(f\) 的 \(n\) 个独立随机变量，并计算出了相应的 \(Y\) 值。通过这种方式，我们生成了 \(r\) 个独立同分布随机变量 \(Y_i = h(X_1^i, X_2^i, \ldots, X_n^i)\)，\(i = 1, \ldots, r\) 的值。
现在，根据大数定律
\[ \lim_{r\rightarrow\infty}\frac{Y_1 + \cdots + Y_r}{r}=E[Y_i]=E[h(X_1, X_2, \ldots, X_n)] \]
所以我们可以使用生成的 \(Y_i\) 值的平均值作为 \(E[h(X_1, X_2, \ldots, X_n)]\) 的估计值。我们称这种近似方法为蒙特卡罗模拟方法。我们每生成一个新的 \(Y\) 值，我们就完成了一次新的模拟运行。
当然，为了使用上述方法，我们需要生成特定密度函数的随机变量。到目前为止，尽管我们只知道如何为均匀随机变量生成特定密度函数的随机变量——使用随机数生成器——但对于自助法（bootstrap）和置换检验（permutation tests）来说已经足够了。因此，接下来的两节将分别介绍自助法和置换检验。然后，在本章的最后几节中，我们将回到如下的问题：如何生成具有任意分布的随机变量以及如何确定何时结束模拟过程。