15  模拟、自助法统计以及置换检验

15.1 引言

在本章中，我们将介绍两种强大的现代统计技术：自助法（bootstrap）统计方法和置换（permutation）检验。这两种方法都是非参数方法，他们不对任何潜在概率分布做特定假设。自助法使我们能够衡量参数估计量的有效性，而置换检验则提供了检验某些统计假设的新方法。然而，这两种方法在实施过程中都需要大量的计算。执行所需计算的最有效且高效的方法是使用 模拟（simulation），这也是本章的第三个主题。

在 章节 15.2 中，我们将介绍模拟的关键——随机数（random numbers），我们将展示如何使用随机数来生成随机排列和随机子集。在 章节 15.2.1 中，我们将介绍用于计算近似期望的蒙特卡罗（Monte Carlo）模拟方法。在 章节 15.3 中，我们将介绍自助法统计，并展示如何应用蒙特卡罗模拟方法来进行所需的分析。在 章节 15.4 中，我们将讨论置换检验：一种用于确定数据序列是否来自单一总体分布的非参数检验方法。在后续的各节中，我们将继续研究模拟。在 章节 15.5 和 章节 15.6 中，我们将展示如何使用随机数来生成任意连续、离散随机变量的值，在 章节 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. 先从数字 \(1,2,\ldots,n\) 中随机选择一个数字，并将该数字放在位置 \(n\)
2. 然后从剩下的 \(n - 1\) 个数字中随机选择一个数字，并将其放在位置 \(n - 1\)
3. 接着从剩下的 \(n - 2\) 个数字中随机选择一个数字，并将其放在位置 \(n - 2\)
4. 以此类推，直到生成整个排列序列（这里的“随机选择”意味着每个可能的选择都有相同的可能性）

然而，我们无需考虑还有哪些确切的数字需要选择，一种方便且有效的方式是：将所有的数字存储在一个有序列表中，然后随机选择有序序列的位置而不是数字本身。也就是说，从数字 \(1,2,\ldots,n\) 的任何排列 \(r_1,r_2,\ldots,r_n\) 开始：

1. 随机选择位置 \(1,\ldots,n\) 中的一个，并将该位置的数字与位置 \(n\) 的数字互换
2. 然后随机选择位置 \(1,\ldots,n - 1\) 中的一个，并将该位置的数字与位置 \(n - 1\) 的数字互换
3. 接着随机选择位置 \(1,\ldots,n - 2\) 中的一个，并将该位置的数字与位置 \(n - 2\) 的数字互换
4. 以此类推

为了实现上述操作，我们需要生成一个随机变量，该随机变量能够等可能地取 \(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\) 的随机排列的算法可以写成如下形式：

1. 令 \(r_1,r_2,\ldots,r_n\) 为数字 \(1,2,\ldots,n\) 的任何一个排列。（例如，我们可以令 \(r_j = j\)，\(j = 1,\ldots,n\)。）
2. 令 \(k = n\)。（要放在位置 \(k\) 的数字有待确定。）
3. 生成一个随机数 \(U\)，并令 \(I=\text{Int}(kU)+1\)。
4. 互换 \(r_I\) 和 \(r_k\) 的值。
5. 令 \(k = k - 1\)。
6. 如果 \(k\gt1\)，转到步骤 3；如果 \(k = 1\)，转到步骤 7。
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\)

Fisher-Yates 洗牌算法

在 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）来说已经足够了。因此，接下来的两节将分别介绍自助法和置换检验。然后，在本章的最后几节中，我们将回到如下的问题：如何生成具有任意分布的随机变量以及如何确定何时结束模拟过程。

15.3 自助法

令 \(X_1, \dots, X_n\) 是来自某个分布（\(F\)）的一个样本，假设我们希望使用该样本来估计该分布的某个参数 \(\theta\)。例如，\(\theta\) 可以是 \(X_i\) 的均值或方差。

假设我们有一个 \(\theta\) 的估计量 \(d = d(X_1, \dots, X_n)\)，我们希望评估 \(d(X_1, \dots, X_n)\) 作为 \(\theta\) 的估计量的优劣性。可以使用均方误差（Mean Square Error, MSE）指标来衡量 \(d(X_1, \dots, X_n)\) 是否可以作为 \(\theta\) 的估计量：

\[ {MSE}_F(d) = E_F[(d(X_1, \dots, X_n) - \theta)^2] \]

也就是说，\(MSE_F(d)\) 是估计量 \(d(X_1, \dots, X_n)\) 与参数 \(\theta\) 之间的距离的平方的期望。这里，我们使用 \(MSE_F\) 和 \(E_F\) 表示：在 \(X_1, \dots, X_n\) 是分布函数为 \(F\) 的独立随机变量的前提下计算的期望。

那么，如何估计呢？

例子 15.3 若 \(\theta = E[X_i]\) 是分布 \(F\) 的均值，且 \(d(X_1,\ldots,X_n)=\overline{X}_n = \sum_{i = 1}^{n}X_i/n\) 是数据值 \(X_1,\ldots,X_n\) 的样本均值，那么因为

\[ E_F[d(X_1,\ldots,X_n)] = E_F[\overline{X}_n]=E_F[X_i]=\theta \]

由此可得

\[ MSE_F(\overline{X}_n)=E_F[(\overline{X}_n - \theta)^2]=\text{Var}_F(\overline{X}_n)=\sigma^2/n \]

其中 \(\sigma^2 = \text{Var}_F(X_i)\)。因此，在这种情况下，可以用样本方差 \(S_n^2\) 来估计总体方差 \(\sigma^2\)，并用 \(S_n^2/n\) 来估计 \(MSE_F(\overline{X}_n)\)，其中

\[ S_n^2=\frac{1}{n - 1}\sum_{i = 1}^{n}(X_i - \overline{X}_n)^2 \]

为 \(X_1,\cdots, X_n\) 的样本方差。 \(\blacksquare\)

在 例子 15.3 中，当使用样本均值作为总体均值的估计量时，我们很容易就能估计其均方误差。然而，如果我们想要估计总体方差呢？也就是说，当 \(\theta = \text{Var}_F(X_i)\) 时将会如何？此时，我们可以使用样本方差来估计总体方差。然而，虽然想到 \(d(X_1,\ldots,X_n)=S_n^2\) 很简单，但如何估计 \(E_F[(d(X_1,\ldots,X_n) - \theta)^2]\) 却并不简单。我们可以使用自助法统计来估计样本反差与总体方差之间的均方误差，接下来我们将对其进行介绍。

为了估计参数 \(\theta\) 与其估计量 \(d(X_1,\ldots,X_n)\) 之间的均方误差，令 \(X_i = x_i\)，\(i = 1,\ldots,n\) 为样本数据值。对于任意 \(x\)，令 \(F_e(x)\) 表示小于或等于 \(x\) 的数据值所占的比例。即，

\[ F_e(x)=\frac{\text{满足 }i\leq n\text{ 且 }x_i\leq x\text{ 的 }i\text{ 的个数}}{n} \]

例如，若 \(n = 5\) ，且 \(X_1 = 5\)，\(X_2 = 3\)，\(X_3 = 9\)，\(X_4 = 2\)，\(X_5 = 6\)，那么

\[ F_e(x)=\begin{cases} 0 & \text{如果 } x < 2 \\ 1/5 & \text{如果 } 2\leq x < 3 \\ 2/5 & \text{如果 } 3\leq x < 5 \\ 3/5 & \text{如果 } 5\leq x < 6 \\ 4/5 & \text{如果 } 6\leq x < 9 \\ 1 & \text{如果 } x\geq 9 \end{cases} \]

我们称函数 \(F_e(x)\) 为经验分布函数（empirical distribution function）。当 \(x_1,\ldots,x_n\) 都不相同时，\(F_e\) 是随机变量 \(X_e\) 的分布函数，\(X_e\) 取 \(x_1,\ldots,x_n\) 中任意一个值的可能性均相等。也就是说，如果数据值都不相同，那么 \(F_e\) 是随机变量 \(X_e\) 的分布函数，使得

\[ P(X_e = x_i)=1/n, \quad i = 1,\ldots,n \]

当数据值存在重复时，\(F_e\) 是随机变量 \(X_e\) 的分布函数，\(X_e\) 取任意指定数据值的概率等于该值在数据集中出现的次数除以 \(n\)。例如，若 \(n = 3\)，且 \(x_1 = x_2 = 1\)，\(x_3 = 2\)，那么 \(X_e\) 是一个随机变量，它取值为 \(1\) 的概率是 \(2/3\)，取值为 \(2\) 的概率是 \(1/3\)。基于对不同值所赋予权重的这种理解，我们仍可以说 \(F_e\) 是一个随机变量的分布函数，该随机变量取 \(x_1, x_2, \ldots, x_n\) 中任意一个值的可能性相等。

现在，对于任意的 \(x\) 值，每个数据值 \(X_i\) （\(i = 1, \ldots, n\) ）小于或等于 \(x\) 的概率为 \(F(x)\)。因此，根据大数定律，当 \(n\) 趋于无穷大时，小于或等于 \(x\) 的数据值的比例以概率 \(1\) 收敛到 \(F(x)\)。所以，当 \(n\) 较大时，\(F_e(x)\) 应接近 \(F(x)\)，这表明经验分布函数 \(F_e\) 可以用作总体分布函数 \(F\) 的一个估计量。

现在，假设 \(\theta_e\) 与分布 \(F_e\) 的关系和 \(\theta\) 与分布 \(F\) 的关系相同。例如，如果 \(\theta\) 是分布函数为 \(F\) 的随机变量 \(X\) 的方差，那么 \(\theta_e\) 就是分布函数为 \(F_e\) 的随机变量 \(X_e\) 的方差。现在，如果 \(F_e\) 接近 \(F\)，那么几乎可以肯定 \(\theta_e\) 会接近 \(\theta\)（从技术层面讲，如果 \(\theta\) 是分布 \(F\) 的连续函数，那么如上的假设就是成立的）。基于这些原因，我们可以近似计算 \(\theta\) 与其估计量 \(d(X_1,\ldots,X_n)\) 之间的均方误差：

\[ MSE_F(d)=E_F[(d(X_1,\ldots,X_n)-\theta)^2]\approx E_{F_e}[(d(X_1,\ldots,X_n)-\theta_e)^2] \]

其中 \(E_{F_e}\) 表示在 \(X_1,\ldots,X_n\) 是相互独立的随机变量，且每个随机变量都具有分布函数 \(F_e\) 的假设下求期望。也就是说，\(X_1,\ldots,X_n\) 中的每一个取 \(x_1,\ldots,x_n\) 中任意一个值的可能性相等。

\[ MSE_{F_e}(d)=E_{F_e}[(d(X_1,\ldots,X_n)-\theta_e)^2] \]

为 \(\theta\) 与其估计量 \(d(X_1,\ldots,X_n)\) 之间的均方误差的自助法估计值（bootstrap estimate of the mean square error）。

现在让我们来看一下，在无需使用 \(MSE(F_e)\) 估计 \(MSE(F)\) 的情况下（当用样本均值来估计总体均值时），\(MSE(F_e)\) 对 \(MSE(F)\) 的估计效果如何。

例子 15.4 对于 例子 15.3，可以使用 \(\overline{X}=\sum_{i=1}^{n}X_i/n\) 作为分布 \(F\) 的均值的估计量。因为随机变量 \(X_e\) 对于每个数据值 \(x_1,\cdots,x_n\) 的权重都是一致的，所以当 \(\theta_e\) 是 \(X_e\) 的均值时，

\[ \theta_e = E[X_e] = \sum_{i=1}^{n}x_iP(X_e=x_i) = \frac{1}{n}\sum_{i=1}^{n}x_i = \overline{x}_n \]

因为

\[ E_{F_e}\left[\sum_{i = 1}^{n}X_i/n\right]=E_{F_e}[X]=\theta_e=\overline{x}_n \]

由此可得

\[ \begin{align*} MSE_{F_e}(\overline{X}_n)&=E_{F_e}\left[\left(\sum_{i = 1}^{n}X_i/n - \overline{x}_n\right)^2\right]\\ &=\text{Var}_{F_e}\left(\sum_{i = 1}^{n}X_i/n\right)\\ &=\frac{1}{n}\text{Var}_{F_e}(X) \end{align*} \]

现在，

\[ \begin{align*} \text{Var}_{F_e}(X)&=E_{F_e}\left[(X - \overline{x}_n)^2\right]\\ &=\sum_{i = 1}^{n}(x_i - \overline{x}_n)^2P_{F_e}(X = x_i)\\ &=\frac{1}{n}\sum_{i = 1}^{n}(x_i - \overline{x}_n)^2 \end{align*} \]

因此，我们已经证明了

\[ MSE_{F_e}(\overline{X}_n)=\frac{1}{n^2}\sum_{i = 1}^{n}(x_i - \overline{x}_n)^2 \]

由于 \(MSE_F(\overline{X}_n)=\frac{1}{n}\text{Var}_F(X)\) 的常用估计量是 \(S_n^2/n\)，其观测值为 \(\frac{1}{n(n - 1)}\sum_{i = 1}^{n}(x_i - \overline{x}_n)^2\)，我们可以看到，自助法估计值与常用估计值几乎相同。\(\blacksquare\)

如前所述，如果数据值为 \(X_i = x_i\)，\(i = 1,\ldots,n\)，那么经验分布函数 \(F_e\) 会给每个点 \(x_i\) 赋予相同的权重 \(1/n\)，因此，通常很容易计算 \(\theta_e\) 的值。为了用自助法来计算 \(\theta\) 和其估计量 \(d(X_1,\ldots,X_n)\) 之间的均方误差，我们需要计算

\[ MSE_{F_e}(d)=E_{F_e}[(d(X_1,\ldots,X_n)-\theta_e)^2] \]

然而，由于如上的期望是在 \(X_1,\ldots,X_n\) 都服从 \(F_e\) 分布的假设下而计算的，所以向量 \((X_1,\ldots,X_n)\) 等可能地取 \(n^n\) 个可能值 \((x_{i_1},x_{i_2},\ldots,x_{i_n})\) 中的任意一个，其中每个 \(i_j\) 是 \(1,\ldots,n\) 中的一个值。因此，除非 \(n\) 很小，否则精确计算 \(MSE_{F_e}(d)\) 是非常困难的。

不过，通过模拟来近似 \(MSE_{F_e}(d)\) 则非常容易。我们生成 \(n\) 个相互独立且服从 \(F_e\) 分布的随机变量 \(X_1^1,\ldots,X_n^1\)，并用它们来计算

\[ Y_1=(d(X_1^1,\ldots,X_n^1)-\theta_e)^2 \]

然后我们重复如上的过程，生成第二组 \(n\) 个相互独立且服从 \(F_e\) 分布的随机变量 \(X_1^2,\ldots,X_n^2\)，并用它们来计算

\[ Y_2=(d(X_1^2,\ldots,X_n^2)-\theta_e)^2 \]

重复该过程很多次，比如 \(r\) 次，得到值 \(Y_1,\ldots,Y_r\)。\(\sum_{i = 1}^{r}Y_i/r\) 将是 \(MSE_{F_e}(d)\) 的近似值，这个近似值随后可作为 \(MSE_F(d)\) 的估计值。

生成一个服从 \(F_e\) 分布的随机变量 \(X\) 是很容易的。只需生成一个随机数 \(U\)，令 \(I = \text{Int}(nU)+1\)，\(I\) 就等可能地取 \(1,\ldots,n\) 中的任意一个值。然后令

\[ X = x_I \]

例子 15.5 假设我们使用来自分布 \(F\) 的样本量为 \(n\) 的样本方差 \(S_n^2=\sum_{i = 1}^{n}(X_i - \overline{X}_n)^2/(n - 1)\) 作为分布 \(F\) 的方差 \(\sigma^2\) 的估计量。为了估计样本方差与 \(\sigma^2\) 的均方误差，设观测数据为 \(X_i = x_i\)，\(i = 1,\ldots,n\)。

因为经验分布函数 \(F_e\) 对所有的值 \(x_i\)（\(i = 1,\ldots,n\)）赋予相等的权重，所以有

\[ E_{F_e}[X]=\sum_{i = 1}^{n}x_iP_{F_e}(X = x_i)=\sum_{i = 1}^{n}x_i/n=\overline{x}_n \]

这表明分布 \(F_e\) 的方差 \(\theta_e\)，可表示为

\[ \theta_e=\text{Var}_{F_e}(X)=E_{F_e}[(X - \overline{x}_n)^2]=\sum_{i = 1}^{n}(x_i - \overline{x}_n)^2/n \]

因此，

\[ MSE_{F_e}(S_n^2)=E_{F_e}[(S_n^2 - \theta_e)^2]=E_{F_e}\left[\left(\frac{\sum_{i = 1}^{n}(X_i - \overline{X}_n)^2}{n - 1}-\theta_e\right)^2\right] \]

我们使用模拟方法来近似 \(MSE_{F_e}(S_n^2)\)。

例如，假设 \(n = 8\)，数据值为 \(x_1 = 5\)，\(x_2 = 9\)，\(x_3 = 12\)，\(x_4 = 8\)，\(x_5 = 7\)，\(x_6 = 15\)，\(x_7 = 3\)，\(x_8 = 6\)。那么

\[ \overline{x}_8 = 8.125, \quad \theta_e=\sum_{i = 1}^{8}(x_i - \overline{x}_8)^2/8\approx13.11 \]

总共会有 \(r\) 次模拟运算，并且 \(N\) 表示当前模拟运算的次数。每次运行中，我们生成 \(8\) 个服从 \(F_e\) 分布的随机变量 \(X_M\)（ \(M = 1,\ldots,8\)）。\(S\) 和 \(SS\) 分别表示到目前为止模拟过程中生成的 \(X_M\) 的总和以及 \(X_M\) 平方的总和。当一次运行完成时，使用以下等式计算样本方差 \(SV\)：

\[ \frac{\sum_{i = 1}^{8}(X_i - \overline{X}_8)^2}{7}=\frac{\sum_{i = 1}^{8}X_i^2 - 8\overline{X}_8^2}{7}=\frac{\sum_{i = 1}^{8}X_i^2 - (\sum_{i = 1}^{8}X_i)^2/8}{7}=\frac{SS - S^2/8}{7} \]

计算 \(SV\) 与 \(\theta_e = 13.11\) 的差值的平方，然后将其加到 \(T\) 中， \(T\) 是之前 \(N - 1\) 个差值平方的总和。当完成 \(r\) 次运算后，模拟结束；样本方差与 \(\theta_e\) 差值平方的平均值就是基于模拟对 \(MSE_{F_e}(S_n^2)\) 的近似值。

1. 令 \(T = 0\)，\(N = 1\)
2. 令 \(M = 1\)
3. \(S = 0\)，\(SS = 0\)
4. 生成一个随机数 \(U\)
5. 令 \(I=\text{Int}(8U)+1\)
6. \(S = S + x_I\)
7. \(SS = SS + x_I^2\)
8. 如果 \(M < 8\) ，令 \(M = M + 1\)，然后转到步骤4
9. \(SV=(SS - S^2/8)/7\)
10. 令 \(T = T+(SV - 13.11)^2\)
11. 如果 \(N < r\)，令 \(N = N + 1\)，然后转到步骤2
12. 如果 \(N = r\) 则返回 \(T/r\) 作为对 \(MSE_{F_e}(S_n^2)\) 的近似值

\(\blacksquare\)

15.4 置换检验

假设我们想要检验原假设 \(H_0\) ：数据 \(X_1, \ldots, X_N\) 是来自某个未指定分布的样本。置换检验 (permutation tests) 是对该假设的一种检验，其中 \(p\) 值的计算是在已知观测到的数据值集合 \(S\) 但不知道具体哪个数据值对应 \(X_1\) 、哪个对应 \(X_2\) 等的前提下进行的。例如，若 \(N = 3\) 且 \(X_1 = 5, X_2 = 7, X_3 = 2\) ，那么 \(p\) 值的计算条件是已知数据值集合 \(S = \{2, 5, 7\}\) 。\(p\) 值的计算利用了这样一个事实：在原假设为真时，给定数据值集合 \(S\) ，将这 \(N\) 个值分配给原始数据的 \(N!\) 种可能方式中的每一种都是等可能的。也就是说，假设 \(N = 3\) 且数据值集合如前所述为 \(S = \{2, 5, 7\}\) 。现在原假设 \(H_0\) 指出 \(X_1, X_2, X_3\) 是独立同分布的。因此，如果 \(H_0\) 为真，那么在给定数据集 \(S\) 的情况下，向量 \((X_1, X_2, X_3)\) 等于值 \(2, 5, 7\) 的 \(3!\) 种置换中的任何一种的概率都是相等的。

置换检验的实施过程如下。根据备择假设，选择一个检验统计量 \(T(X_1, \ldots, X_N)\) 。暂时假设检验统计量的较大值是备择假设的证据。然后观测数据值，设 \(X_i = x_i, i = 1, \ldots, N\) ，并计算 \(T(x_1, \ldots, x_N)\) 的值。令 \(S = \{x_1, \ldots, x_N\}\) 为由 \(N\) 个观测值组成的无序集合。那么，如果检验统计量的值为 \(T(x_1, \ldots, x_N) = t\) ，则由这些数据产生的原假设的 \(p\) 值为：

\[ p\text{-value} = P_{H_0}(T(X_1, \ldots, X_N) \geq t \mid S = \{x_1, \ldots, x_N\}) \tag{15.1}\]

现在，在 \(H_0\) 下，\(X_1, \ldots, X_N\) 等可能地等于 \(x_1, \ldots, x_N\) 的 \(N!\) 种置换中的任何一种。因此，令 \(I_1, \ldots, I_N\) 为一个随机向量，它等可能地取 \(1, \ldots, N\) 的 \(N!\) 种置换中的任何一种，我们可以将上述 \(p\) 值写为：

\[ \begin{align*} p\text{-value} &= P(T(x_{I_1}, x_{I_2}, \ldots, x_{I_N}) \geq t) \\ &= \frac{\text{满足 } T(x_{i_1}, x_{i_2}, \ldots, x_{i_N}) \geq t \text{ 的置换 } (i_1, \ldots, i_N) \text{ 的个数}}{N!} \end{align*} \tag{15.2}\]

为了说明，假设要在 \(N\) 周内观测数据，其中 \(X_i\) 是在第 \(i\) 周观测到的数据值，并且我们想要检验原假设：

\(H_0\) ：\(X_1, \ldots, X_N\) 是独立同分布的

对抗备择假设：

\(H_1\) ：\(X_i\) 随 \(i\) 的增加而呈现增加趋势

现在，如果原假设为真且数据是独立同分布的，那么在已知值集合 \(\{X_1, \ldots, X_N\}\) 但不知道哪个值对应 \(X_1\) 或哪个对应 \(X_2\) 等的情况下，统计量 \(\sum_{j=1}^N j X_j\) 的分布就好像我们随机配对了两个数据集 \(\{1, \ldots, N\}\) 和 \(\{X_1, \ldots, X_N\}\) 然后对这 \(N\) 对配对值的乘积求和。另一方面，如果备择假设为真，那么 \(\sum_{j=1}^N j X_j\) 往往比我们仅仅随机配对值时的趋势更大。这是因为两个同等大小集合的配对值之和在最大值与最大值配对、次大值与次大值配对等情况下达到最大。（在统计学术语中，当 \(X_j\) 随 \(j\) 增加而增加时，数据对 \((j, X_j), j = 1, \ldots, N\) 的相关系数较大。）因此，\(H_0\) 对抗 \(H_1\) 的一种可能的置换检验是：

1. 观测数据值——设 \(X_j = x_j, j = 1, \ldots, N\) 。
2. 令 \(t = \sum_{j=1}^N j x_j\) 。
3. 确定由下式给出的 \(p\) 值： \[ p\text{-value} = P\left(\sum_{j=1}^N I_j x_j \geq t\right) \tag{15.3}\] 其中 \(I_1, \ldots, I_N\) 等可能地为 \(1, \ldots, N\) 的 \(N!\) 种置换中的任何一种。

上述 \(p\) 值可以通过使用 例子 15.2 的方法生成随机置换进行模拟近似。

例子 15.6 为了确定 DVD 播放机的周销量是否呈下降趋势，一家大型电子商店的经理记录了过去 12 周的此类销量，从第 1 周到第 12 周（本周）的结果如下： \[22, 24, 20, 18, 16, 14, 15, 15, 13, 17, 12, 14\] 这些数据是否足够强大，足以拒绝销量分布随时间不变的原假设，从而使经理能够得出销量呈下降趋势的结论？

解 ：令原假设为销量分布随时间不变，备择假设为销量呈下降趋势。因此，如果备择假设为真，那么第 \(j\) 周的销量 \(X_j\) 与 \(j\)之间将存在负相关。因此，\(\sum_{j=1}^{12} j X_j\) 的相对较小值将是支持备择假设的证据。现在，设 \(x_j\) 为 \(X_j\) 的观测值，销售数据给出： \[ \sum_{j=1}^{12} j x_j = 1178 \] 因此，数据来自同一分布的原假设对抗数据随时间递减的备择假设的置换检验的 \(p\) 值由下式给出： \[ p\text{-value} = P\left(\sum_{j=1}^{12} I_j x_j \leq 1178\right) \] 其中 \(I_1, \ldots, I_{12}\) 等可能地为 \(1, \ldots, 12\) 的 \(12!\) 种置换中的任何一种。使用 \(10^5\) 次运行进行的模拟得出： \[ p\text{-value} \approx .00039 \] 这引导我们拒绝销量分布随时间不变的原假设。 \(\blacksquare\)

虽然 \(\sum_{j=1}^N j X_j\) 是最常用于检验 \(X_1, \ldots, X_n\) 独立同分布原假设对抗 \(X_j\) 随 \(j\) 增加而增加的备择假设的检验统计量，但它并不是唯一的可能性。实际上，我们可以选择任何形式为 \(\sum_{j=1}^N a_j X_j\) 的检验统计量，其中 \(a_1 < a_2 < \ldots < a_n\) 。（例如，我们可以选择 \(a_j = j^2\) 。）与前述类似，首先确定统计量的值，设其为 \(t\) 。因为当备择假设为真时，较大的 \(a_j\) 值往往会与较大的数据值配对，从而使得 \(\sum_{j=1}^N a_j X_j\) 趋向于变得更大。因此，当 \(t\) 较大时，我们同样希望拒绝原假设。于是，所得的 \(p\) 值为： \[ p\text{-value} = P\left(\sum_{j=1}^N a_{I_j} x_j \geq t\right) \tag{15.4}\] 其中 \(I_1, \ldots, I_N\) 等可能地为 \(1, \ldots, N\) 的 \(N!\) 种置换中的任何一种。

根据备择假设的不同，我们可以选择其他的常数 \(a_j, j = 1, \ldots, N\) 来检验数据值独立同分布的原假设。例如，如果备择假设是数据在中间值趋于更高而在极端值趋于更低，那么我们可以让检验统计量的形式为 \(T = \sum_{j=1}^N a_j X_j\) ，其中 \(a_1, \ldots, a_N\) 使得中间位置的常数比两端的更大。例如，我们可以使用 \(a_j = j(N - j), j = 1, \ldots, N\) 。由于这同样使得当备择假设为真时，较大的数据值更有可能与较大的常数配对，因此当 \(T\) 较大时，我们仍然希望拒绝原假设。

15.4.1 置换检验中的正态近似

虽然不如进行模拟那样准确，但可以通过假设检验统计量近似服从正态分布来估算置换检验的 \(p\) 值。现在，在数据值独立同分布的原假设下，给定数据集 \(S = \{x_1, \ldots, x_N\}\) ，随机变量 \(X_i\) 等可能地为这 \(N\) 个值中的任何一个，并且当 \(i \neq j\) 时，随机向量 \((X_i, X_j)\) 等可能地取 \(x_k, x_r (r \neq k)\) 这 \(N(N - 1)\) 个值中的任何一个。因此，给定 \(S = \{x_1, \ldots, x_N\}\) ：

\[ E[X_i] = \frac{1}{N} \sum_{i=1}^N x_i = \bar{x} \tag{15.5}\] \[ E[X_i^2] = \frac{1}{N} \sum_{i=1}^N x_i^2 \tag{15.6}\] \[ \begin{align*} E[X_i X_j] &= \frac{1}{N(N - 1)} \sum_k \sum_{r \neq k} x_k x_r \\ &= \frac{1}{N(N - 1)} \left( \sum_k \sum_r x_k x_r - \sum_k \sum_{r=k} x_k x_r \right) \\ &= \frac{1}{N(N - 1)} \left( \sum_k x_k \sum_r x_r - \sum_k x_k^2 \right) \\ &= \frac{1}{N(N - 1)} \left( N^2 \bar{x}^2 - \sum_{k=1}^N x_k^2 \right) \end{align*} \tag{15.7}\]

因此，令 \(v = \text{Var}(X_i)\) 且 \(c = \text{Cov}(X_i, X_j), i \neq j\) ，上述推导给出：

\[ E[X_i] = \bar{x} \] \[ v = \text{Var}(X_i) = \frac{1}{N} \sum_{i=1}^N x_i^2 - \bar{x}^2 \tag{15.8}\] \[ c = \text{Cov}(X_i, X_j) = \frac{1}{N(N - 1)} (N^2 \bar{x}^2 - \sum_{k=1}^N x_k^2) - \bar{x}^2 \tag{15.9}\] \[ v - c = \frac{\sum_{i=1}^N x_i^2 - N \bar{x}^2}{N - 1} \tag{15.10}\]

因此，当 \(H_0\) 为真时，检验统计量 \(T = \sum_{j=1}^N j X_j\) 的均值为： \[ E[T] = \frac{N(N + 1)}{2} \bar{x} \tag{15.11}\] 方差为： \[ \begin{align*} \text{Var}(T) &= \text{Var}\left(\sum_{j=1}^N j X_j\right) \\ &= \sum_{j=1}^N \text{Var}(j X_j) + \sum_i \sum_{j \neq i} \text{Cov}(i X_i, j X_j) \\ &= v \sum_{j=1}^N j^2 + c \sum_i \sum_{j \neq i} i j \\ &= v \sum_{j=1}^N j^2 + c \left( \sum_i \sum_j i j - \sum_i \sum_{j=i} i j \right) \\ &= v \sum_{j=1}^N j^2 + c \left( \sum_{i=1}^N i \sum_{j=1}^N j - \sum_{i=1}^N i^2 \right) \\ &= (v - c) \sum_{j=1}^N j^2 + \frac{c N^2(N + 1)^2}{4} \\ &= (v - c) \frac{N(N + 1)(2N + 1)}{6} + \frac{c N^2(N + 1)^2}{4} \end{align*} \tag{15.12}\]

利用上述结果，我们可以通过假设 \(T\) 在 \(H_0\) 为真时近似服从正态分布，来估算置换检验的 \(p\) 值。

例子 15.7 再次考虑 例子 15.6 。计算得出，在 \(H_0\) 下： \[ E[T] = 1300 \quad \text{Var}(T) = 1958.81 \] 因此，正态近似得出： \[ \begin{align*} p\text{-value} &= P_{H_0}(T \leq 1178) \\ &= P_{H_0}\left( \frac{T - 1300}{\sqrt{1958.81}} \leq \frac{1178 - 1300}{\sqrt{1958.81}} \right) \\ &\approx \Phi(-2.757) = .0029 \end{align*} \] 这与模拟给出的值非常接近。

现在假设 12 个数据值虽然与之前相同，但它们现在的出现顺序为： \[22, 14, 14, 16, 24, 20, 18, 15, 17, 15, 12, 13\] 对于这些数据，检验统计量的值为 \(\sum_{j=1}^{12} j X_j = 1233\) ，正态近似得出： \[ \begin{align*} p\text{-value} &= P_{H_0}(T \leq 1233) \\ &= P_{H_0}\left( \frac{T - 1300}{\sqrt{1958.81}} \leq \frac{1233 - 1300}{\sqrt{1958.81}} \right) \\ &\approx \Phi(-1.514) = .065 \end{align*} \]

最后，再次假设这 12 个数据值与之前相同，但假设它们现在的出现顺序为： \[22, 14, 14, 16, 24, 13, 18, 15, 17, 15, 12, 20\] 在这种情况下，检验统计量的值为 \(\sum_{j=1}^{12} j X_j = 1275\) 。因此，正态近似得出： \[ \begin{align*} p\text{-value} &= P_{H_0}(T \leq 1275) \\ &= P_{H_0}\left( \frac{T - 1300}{\sqrt{1958.81}} \leq \frac{1275 - 1300}{\sqrt{1958.81}} \right) \\ &\approx \Phi(-.565) \approx .286 \end{align*} \] \(10^5\) 次运行的模拟得出了与上述非常相似的值。模拟给出： \[ P_{H_0}(T \leq 1233) \approx .068 \] 且 \[ P_{H_0}(T \leq 1275) \approx .299 \] 这些值与正态近似给出的值相当接近。 \(\blacksquare\)

例子 15.8 为了进一步说明正态近似的有效性，假设 \(N = 4\) ，且数据按以下顺序出现： \[13, 7, 5, 3\] 假设我们想使用这些数据来检验原假设：数据是来自某个分布的样本，对抗备择假设：数据趋于递减。检验统计量的值为 \(T = \sum_{j=1}^4 j X_j = 54\) 。简单的计算得出： \[ c = -4.667, \quad v = 14 \] 表明在 \(H_0\) 下： \[ E[T] = 70, \quad \text{Var}(T) = 93.33 \] 因此，令 \(Z\) 为标准正态随机变量，正态近似得出： \[ \begin{align*} p\text{-value} &= P_{H_0}(T \leq 54) \\ &= P_{H_0}\left( \frac{T - 70}{\sqrt{93.33}} \leq \frac{54 - 70}{\sqrt{93.33}} \right) \\ &\approx P(Z \leq -1.656) = .049 \end{align*} \] 而精确值为： \[ p\text{-value} = P_{H_0}(T \leq 54) = 1/4! \approx .042 \] \(\blacksquare\)

15.4.2 两样本置换检验

置换检验在两样本问题中也很有用，即检验来自两个总体的样本是否具有相同的潜在分布。具体而言，令 \(X_1, \ldots, X_n\) 是来自未知总体分布 \(F\) 的样本，并令 \(X_{n+1}, \ldots, X_{n+m}\) 是来自未知总体分布 \(G\) 的独立样本。假设我们要利用这些数据来检验原假设：两个总体分布是相同的，对抗备择假设：来自第二个分布的数据往往比来自第一个分布的数据更大。也就是说，我们要检验原假设：

\(H_0: F = G\)

对抗备择假设：

\(H_1\) ：来自 \(G\) 的数据倾向于比来自 \(F\) 的数据更大

如果数据值为 \(X_i = x_i, i = 1, \ldots, n + m\) ，那么上述原假设的置换检验是在已知 \(S = \{x_1, \ldots, x_{n+m}\}\) （这 \(n+m\) 个数且无特定顺序）的前提下进行的。如果 \(H_0\) 为真，那么 \(n+m\) 个随机变量 \(X_1, \ldots, X_{n+m}\) 都是独立同分布的，那么给定值集合 \(S\) ，该集合中大小为 \(n\) 的每个子集作为 \(X_1, \ldots, X_n\) 的数据值集合的可能性是相等的。由于备择假设是总体分布 \(F\) 的数据趋向于小于总体分布 \(G\) 的数据，一个合理的检验是：如果从数据集 \(S\) 中随机选择 \(n\) 个值时，来自总体分布 \(F\) 的数据值之和小于可能由机会预期的值，则拒绝原假设。更具体地说，我们可以通过计算 \(\sum_{i=1}^n x_i\) 来检验 \(H_0\) ；设其值为 \(t\) 。那么 \(H_0\) 对抗 \(H_1\) 的置换检验的 \(p\) 值将等于从值 \(x_1, \ldots, x_{n+m}\) 中随机选择 \(n\) 个值的和小于或等于 \(t\) 的概率。即：

\[ p\text{-value} = P\left( \sum_{i \in R} x_i \leq t \right) \tag{15.13}\]

其中 \(R\) 等可能地为从集合 \(\{1, 2, \ldots, n + m\}\) 中选出的 \(\binom{n+m}{n}\) 个大小为 \(n\) 的子集中的任何一个。虽然仅当 \(\binom{n+m}{n}\) 较小时才可能进行精确计算，但通过模拟很容易获得精确的近似值。在每次模拟运行中，我们使用 例子 15.1 的方法随机生成值 \(1, \ldots, n + m\) 的一个大小为 \(n\) 的子集。如果得到的子集是 \(R\) ，那么我们检查 \(\sum_{i \in R} x_i\) 是否小于或等于 \(t\) 。模拟运行中满足该条件的比例即为我们对前述 \(p\) 值的估计。

注 ：我们再次可以使用正态近似，而不是模拟，来估计 \(p\) 值。从随机子集 \(R\) 开始（等可能地为集合 \(\{1, 2, \ldots, n + m\}\) 的 \(\binom{n+m}{n}\) 个大小为 \(n\) 的子集中的任何一个），令对于 \(i = 1, \ldots, n + m\) ：

\[ I_i = \begin{cases} 1, & \text{如果 } i \in R \\ 0, & \text{如果 } i \notin R \end{cases} \]

那么， \[ \sum_{i \in R} x_i = \sum_{i=1}^{n+m} x_i I_i \tag{15.14}\]

通过与 章节 15.4.1 中类似的分析，我们现在可以得出：

\[ E_{H_0}\left[ \sum_{i \in R} x_i \right] = E_{H_0}\left[ \sum_{i=1}^{n+m} x_i I_i \right] = n \bar{x} \tag{15.15}\]

且

\[ \text{Var}_{H_0}\left( \sum_{i \in R} x_i \right) = \text{Var}_{H_0}\left( \sum_{i=1}^{n+m} x_i I_i \right) = \frac{nm}{n + m - 1} \left( \frac{\sum_{i=1}^{n+m} x_i^2}{n + m} - \bar{x}^2 \right) \tag{15.16}\]

其中 \(\bar{x} = \sum_{i=1}^{n+m} x_i / (n + m)\) 。

15.5 离散随机变量的生成

假设我们想要生成一个具有如下概率质量函数的随机变量 \(X\) 的值： \[ P(X = x_i) = p_i, \quad i = 1, \ldots, \quad \sum_i p_i = 1 \tag{15.17}\]

为了生成 \(X\) 的值，我们先生成一个随机数 \(U\) ，并设定： \[ X = x_i \quad \text{如果 } p_1 + \ldots + p_{i-1} < U \leq p_1 + \ldots + p_{i-1} + p_i \tag{15.18}\]

也就是说， \[ X = \begin{cases} x_1, & \text{如果 } U \leq p_1 \\ x_2, & \text{如果 } p_1 < U \leq p_1 + p_2 \\ x_3, & \text{如果 } p_1 + p_2 < U \leq p_1 + p_2 + p_3 \\ \vdots \\ x_i, & \text{如果 } p_1 + \ldots + p_{i-1} < U \leq p_1 + \ldots + p_i \\ \vdots \end{cases} \tag{15.19}\]

由于 \(U\) 在 \((0, 1)\) 上均匀分布，因此对于 \(0 < a < b < 1\) ，有 \[ P(a < U \leq b) = b - a \tag{15.20}\] 从而有， \[ P\left( \sum_{j=1}^{i-1} p_j < U \leq \sum_{j=1}^i p_j \right) = p_i \tag{15.21}\] 这表明 \(X\) 具有所需的概率质量函数。这种生成 \(X\) 的方法被称为离散逆变换方法 (discrete inverse transform method)。

例子 15.9 要生成一个伯努利随机变量 \(X\) ，使得 \[ P(X = 1) = p = 1 - P(X = 0) \] 生成一个随机数 \(U\) ，并设定 \[ X = \begin{cases} 1, & \text{如果 } U \leq p \\ 0, & \text{如果 } U > p \end{cases} \] \(\blacksquare\)

例子 15.10 假设现在我们想要生成一个参数为 \(n\) 和 \(p\) 的二项随机变量 \(X\) 。回想一下，\(X\) 表示 \(n\) 次独立试验中成功的次数，每次试验成功的概率为 \(p\) 。我们可以通过生成这 \(n\) 次试验的结果来生成 \(X\) 。也就是说，我们可以生成 \(n\) 个随机数 \(U_1, \ldots, U_n\) ，如果 \(U_i \leq p\) ，则称试验 \(i\) 成功，然后设定 \[ X = \text{满足 } U_i \leq p \text{ 的 } i \text{ 的个数} \] 另一种可能性是使用逆变换方法。

为了有效地使用逆变换方法，我们需要一种有效的方法来递归地计算以下值： \[ p_i = P(X = i) = \binom{n}{i} p^i (1 - p)^{n-i}, \quad i = 0, \ldots, n \tag{15.22}\] 这可以通过首先注意到以下等式来实现： \[ \frac{\binom{n}{i+1}}{\binom{n}{i}} = \frac{n! / [(n - i - 1)! (i + 1)!]}{n! / [(n - i)! i!]} = \frac{n - i}{i + 1} \tag{15.23}\] 从而得出： \[ \frac{p_{i+1}}{p_i} = \frac{n - i}{i + 1} \frac{p^{i+1} (1 - p)^{n-i-1}}{p^i (1 - p)^{n-i}} = \frac{n - i}{i + 1} \frac{p}{1 - p} \tag{15.24}\] 因此， \[ p_{i+1} = \frac{n - i}{i + 1} \frac{p}{1 - p} p_i \tag{15.25}\]

利用上述结果，我们现在可以给出生成二项 \((n, p)\) 随机变量 \(X\) 的逆变换方法。在下文中，\(i\) 表示 \(X\) 的可能值，变量 \(P\) 是 \(X = i\) 的概率，而变量 \(F\) 是 \(X \leq i\) 的概率。（也就是说，对于给定的 \(i\)，\(P = p_i\) 且 \(F = \sum_{j=0}^i p_j\) 。）此外，令 \(\alpha = p_0 = (1 - p)^n\) ，并令 \(b = p / (1 - p)\) 。

1. 设定 \(i = 0, P = \alpha, F = \alpha\)
2. 生成一个随机数 \(U\)
3. 如果 \(U \leq F\) ，设定 \(X = i\) 并停止
4. \(P = \frac{n - i}{i + 1} b P\)
5. \(F = F + P\)
6. \(i = i + 1\)
7. 转到第 3 步

（在前述过程中，当我们要执行 \(P = \frac{n - i}{i + 1} b P\) 时，并不是从字面上理解为一个代数等式；而是指将 \(P\) 的值更新。它的新值是其旧值乘以 \(\frac{n - i}{i + 1} b\) 。类似地，当我们写 \(F = F + P\) 时，是指通过将 \(P\) 加到 \(F\) 的旧值上来更新 \(F\) 的值。）

由于该算法首先检查是否 \(X = 0\) ，然后检查是否 \(X = 1\) ，依此类推，因此所需的迭代次数（即进入第 3 步的次数）比 \(X\) 的最终值多一次。因此，平均而言，该算法需要 \(E[X + 1] = np + 1\) 次迭代来生成 \(X\) 的值。 \(\blacksquare\)

15.6 连续随机变量的生成

令 \(F\) 为连续随机变量的分布函数。对于介于 \(0\) 和 \(1\) 之间的任何 \(u\) ，量 \(F^{-1}(u)\) 被定义为满足 \(F(x) = u\) 的 \(x\) 值。也就是说，\(F(F^{-1}(u)) = u\) 。由于连续随机变量的分布函数是严格递增的，因此 \(F^{-1}(u)\) 的值是唯一的。我们称 \(F^{-1}\) 为 \(F\) 的反函数。

一种生成具有分布函数 \(F\) 的连续随机变量的通用方法被称为逆变换方法 (inverse transformation method)，它基于以下命题。

定理 15.1 令 \(U\) 为均匀 \((0, 1)\) 随机变量。对于任何连续分布函数 \(F\) ，如果我们定义 \[ X = F^{-1}(U) \tag{15.26}\] 那么 \(X\) 具有分布函数 \(F\) 。

证明 ：因为分布函数 \(F\) 是非减的，所以对于任何数字 \(a\) 和 \(b\) ，不等式 \(a \leq b\) 等价于不等式 \(F(a) \leq F(b)\) 。因此， \[ P(F^{-1}(U) \leq x) = P(F(F^{-1}(U)) \leq F(x)) = P(U \leq F(x)) = F(x) \tag{15.27}\] 从而表明 \(F^{-1}(U)\) 具有分布 \(F\) 。 \(\blacksquare\)

例子 15.11 （生成指数随机变量）令 \[ F(x) = 1 - e^{-\lambda x}, \quad x \geq 0 \tag{15.28}\] 为具有参数 \(\lambda\) 的指数随机变量的分布函数。那么 \(F^{-1}(u)\) 是满足下式的 \(x\) 值： \[ u = F(x) = 1 - e^{-\lambda x} \] 或等价地， \[ e^{-\lambda x} = 1 - u \] 或 \[ -\lambda x = \log(1 - u) \] 或 \[ x = -\frac{1}{\lambda} \log(1 - u) \] 因此，根据 定理 15.1，我们可以通过生成一个均匀 \((0, 1)\) 随机变量 \(U\) 并设定如下方式来生成具有参数 \(\lambda\) 的指数随机变量 \(X\) ： \[ X = -\frac{1}{\lambda} \log(1 - U) \tag{15.29}\] 由于 \(1 - U\) 也是一个均匀 \((0, 1)\) 随机变量，因此得出 \(-\frac{1}{\lambda} \log(1 - U)\) 和 \(-\frac{1}{\lambda} \log(U)\) 具有相同的分布，从而表明 \[ X = -\frac{1}{\lambda} \log(U) \tag{15.30}\] 也是参数为 \(\lambda\) 的指数分布。 \(\blacksquare\)

15.6.1 正态随机变量的生成

由于反转正态随机变量的分布函数在计算上比较复杂，因此使用了特殊的方法来生成正态随机变量。以下方法被称为 Box-Muller 方法 (Box-Muller method)。

首先，假设 \(X\) 和 \(Y\) 是独立的标准正态随机变量，因此它们的联合密度函数为 \[ f(x, y) = \frac{1}{\sqrt{2\pi}} e^{-x^2 / 2} \frac{1}{\sqrt{2\pi}} e^{-y^2 / 2} = \frac{1}{2\pi} e^{-(x^2 + y^2) / 2}, \quad -\infty < x, y < \infty \tag{15.31}\] 令 \(R, \Theta\) 为点 \((X, Y)\) 的极坐标 (polar coordinates)。现在根据定义，\(R^2 = X^2 + Y^2\) 是具有 2 个自由度的卡方随机变量 (chi-square random variable)，如 章节 5.8.1 所示，该分布与参数为 \(1/2\) 的指数分布 (exponential distribution) 相同（即均值为 2）。因此，\(R^2\) 的密度函数为 \[ f_{R^2}(r) = \frac{1}{2} e^{-r / 2}, \quad 0 < r < \infty \tag{15.32}\] 现在考虑给定 \(R^2 = r\) 时 \(X, Y\) 的条件联合密度函数。因为 \[ f(x, y) = \frac{1}{2\pi} e^{-r / 2} \quad \text{当 } x^2 + y^2 = r \] 是一个常数（当 \(x^2 + y^2 = r\) 时），因此直观上可以证明（并且可以证明）：在给定 \(R^2 = r\) 的条件下，向量 \((X, Y)\) 均匀分布在半径为 \(\sqrt{r}\) 的圆周上。但这暗示了：在给定 \(R^2 = r\) 的条件下，点 \((X, Y)\) 的极坐标 \(\Theta\) 在 \((0, 2\pi)\) 上均匀分布。因为这对所有的 \(r\) 都成立，所以极坐标 \(R\) 和 \(\Theta\) 是独立的，其中 \(R\) 服从均值为 2 的指数随机变量的平方根分布，而 \(\Theta\) 是 \((0, 2\pi)\) 上的均匀随机变量。

利用前述内容，我们可以通过首先生成它们的极坐标 \(R\) 和 \(\Theta\) 来生成独立的标准正态随机变量 \(X\) 和 \(Y\) 。由于 \(-\log(U)\) 是均值为 1 的指数分布，我们可以通过生成独立的均匀 \((0, 1)\) 随机变量 \(U_1\) 和 \(U_2\) ，然后设定如下方式来生成 \((X, Y)\) 的极坐标： \[ R^2 = -2 \log(U_1) \] 且 \[ \Theta = 2\pi U_2 \] 利用极坐标 \(R, \Theta\) 回到直角坐标的公式 \[ X = R \cos(\Theta), \quad Y = R \sin(\Theta) \] 表明 \[ X = \sqrt{-2 \log(U_1)} \cos(2\pi U_2) \tag{15.33}\] \[ Y = \sqrt{-2 \log(U_1)} \sin(2\pi U_2) \tag{15.34}\] 是独立的标准正态随机变量。

要生成均值为 \(\mu\) 、方差为 \(\sigma^2\) 的正态随机变量，只需生成独立的标准正态变量 \(X\) 和 \(Y\) ，然后取变量 \(\mu + \sigma X\) 和 \(\mu + \sigma Y\) 。

15.7 确定蒙特卡罗研究中的模拟运行次数

假设我们要生成 \(r\) 个独立同分布的随机变量 \(Y_1, \ldots, Y_r\) ，它们的均值为 \(\mu\) ，以便使用 \[ \bar{Y}_r = \sum_{i=1}^r Y_i / r \tag{15.35}\] 作为 \(\mu\) 的估计量。现在，设 \(\sigma^2\) 为 \(Y_i\) 的方差，由中心极限定理可知，\(\bar{Y}_r\) 将近似服从均值为 \(\mu\) 、方差为 \(\sigma^2/r\) 的正态分布。因此，我们可以 95% 肯定 \(\mu\) 落在区间 \[ (\bar{Y}_r - 1.96 \sigma / \sqrt{r}, \quad \bar{Y}_r + 1.96 \sigma / \sqrt{r}) \tag{15.36}\] 内。（更一般地，我们可以 \(100(1 - \alpha)\%\) 肯定 \(\mu\) 将在 \(\bar{Y}_r \pm z_{\alpha/2} \sigma / \sqrt{r}\) 之间。）

因此，如果 \(\sigma^2\) 已知，我们可以选择 \(r\) 来达到我们所需的准确度水平。然而，几乎总是会出现这样的情况：\(\sigma^2\) 和 \(\mu\) 一样，都是未知的。为了克服这一困难，我们可以进行两阶段模拟实验。在第一阶段，我们生成 \(k\) 次运行，其中 \(k\) 通常远小于研究中预期使用的次数。进行这些运行会产生随机变量 \(Y_1, \ldots, Y_k\) 的值。然后我们使用这些值的样本方差 \[ S_k^2 = \frac{1}{k - 1} \sum_{i=1}^k (Y_i - \bar{Y}_k)^2 \tag{15.37}\] 来估计 \(\sigma^2\) 。然后，假设这就是 \(\sigma^2\) 的真实值，我们确定 \(r\) 的合适值。接着，在模拟的第二阶段，我们再生成额外的 \(r - k\) 次运行。

Problems