50、随机采样技术全解析

随机采样技术全解析

1. 随机采样概述

许多算法都会用到随机数，这就要求我们能根据特定概率密度 $p(x)$ 从集合中选取元素 $x$。多次重复选取后，特定元素 $\tilde{x}$ 出现的频率应与概率 $p(\tilde{x})$ 成正比。下面将介绍从连续和离散随机变量中采样的通用技术。

2. 随机数生成器

2.1 真正随机数与伪随机数

计算机一般无法生成真正的随机数，原因有二：一是数字计算机只能用有限位数近似表示实数；二是计算机算法是确定性的，只能产生有限且可预测的输出，而真正的随机序列是无限且不可重现的。

在实际应用中，伪随机数序列通常就足够了。伪随机数由确定性程序生成，但具有随机数的一些关键特征：
- 序列值在 $[0, 1]$ 上均匀分布。
- 序列元素不相关。
- 即使知道已生成的所有元素，也很难猜出下一个元素的值。

2.2 伪随机数序列的验证

验证伪随机数序列 ${\xi_1, \xi_2, \ldots, \xi_R}$ 在 $[0, 1]$ 上的均匀分布相对容易，只需生成大量元素 $R$ 并绘制频率直方图，随着 $R$ 增大，直方图应快速收敛到均匀分布。检查序列中是否存在相关性稍难，原则上要检查任意阶的相关性，但实际上，两两不相关意味着 $(\xi_n, \xi_{n + 1})$ 在 $[0, 1]^2$ 上均匀分布，三项不相关意味着 $(\xi_n, \xi_{n + 1}, \xi_{n + 2})$ 在单位立方体上均匀分布。

验证伪随机数生成器生成的序列是否不可预测相对困难，大多数经典随机数生成器在特定参数值下会产生可预测的序列。

2.3 线性同余生成器（LCG）

线性同余生成器（LCG）是一种经典的随机数生成器，由递归关系定义：
$x_{n + 1} = (ax_n + c) \bmod m, n > 0$
其中 $x_0, a, c, m \in N$。该方程生成 $0$ 到 $m - 1$ 之间的整数序列 $x_1, x_2, \ldots$，将 $x_n$ 除以 $m$ 得到有理数 $\xi_n = x_n/m \in [0, 1 - 1/m]$。

若参数 $x_0, a, c, m$ 选择得当，计算得到的序列 $\xi_n$ 近似于 $[0, 1]$ 上的均匀随机变量。$\xi_n$ 的周期至多为 $m$，因此 $m$ 要足够大，以生成周期长于所需随机数数量的序列。一般来说，序列质量取决于 $a, c, m$ 的选择，当且仅当满足以下三个条件时，LCG 生成的序列周期恰好为 $m$：
- $m$ 和 $c$ 互质，即 $m$ 和 $c$ 没有公因数。
- $(a - 1)$ 能被 $m$ 的所有质因数整除。
- 若 $m$ 是 $4$ 的倍数，则 $(a - 1)$ 是 $4$ 的倍数。

两个较好的参数选择是 $m = 2^{32}, a = 1,664,525, c = 1,013,904,223$ 或 $m = 2^{31} - 1, a = 75, c = 0$。需要注意的是，虽然线性同余生成器实现简单，但并非生成伪随机数序列的最佳方法。目前，梅森旋转算法（Mersenne Twister）是大多数编程语言和科学计算库的首选伪随机数生成器，其过程较为复杂。

2.4 流程图：线性同余生成器工作流程

graph TD;
    A[初始化 x0, a, c, m] --> B[计算 x1 = (ax0 + c) mod m];
    B --> C[计算 ξ1 = x1/m];
    C --> D{是否需要更多随机数};
    D -- 是 --> E[计算 xn+1 = (axn + c) mod m];
    E --> F[计算 ξn+1 = xn+1/m];
    F --> D;
    D -- 否 --> G[结束];

3. 逆函数法

3.1 原理

设 $x \in R$ 是具有概率密度函数 $p_x(x)$ 的随机变量，即 $p_x(x) \geq 0 \forall x \in R$ 且 $\int_{R} p_x(x) dx = 1$。若要根据 $p_x(x)$ 采样元素，逆函数法是一种简单方法。该方法基于累积分布函数的定义：
$P_x(x) = Pr(x \leq x) = \int_{-\infty}^{x} p_x(t) dt$
设随机变量 $y = P_x(x)$ 且在 $[0, 1]$ 上均匀分布。若在 $[0, 1]$ 中随机采样一个值 $u$，则 $\xi = P_x^{-1}(u)$ 按 $p_x(x)$ 分布。

3.2 指数分布采样

若要使用逆函数法从指数分布 $p_x(x) = \lambda e^{-\lambda x}, x \geq 0$ 中采样，累积分布函数为：
$P_x(x) = \lambda \int_{0}^{x} e^{-\lambda t} dt = 1 - e^{-\lambda x}$
其逆函数为：
$P_x^{-1}(u) = -\frac{1}{\lambda} \log (1 - u)$
因此，在 $[0, 1]$ 上均匀采样 $u$，通过计算 $\xi = -\frac{1}{\lambda} \log (1 - u)$ 可得到按 $p_x(x) = \lambda e^{-\lambda x}$ 分布的样本。由于 $u$ 和 $1 - u$ 都在 $[0, 1]$ 上均匀分布，也可计算 $\xi = -\frac{1}{\lambda} \log(u)$。

3.3 连续幂律分布采样

若 $p_x(x)$ 是幂律分布 $p_x(x) = (\gamma - 1)x^{-\gamma}, x \geq 1, \gamma > 1$，累积分布函数为：
$P_x(x) = (\gamma - 1) \int_{1}^{x} t^{-\gamma} dt = 1 - x^{1 - \gamma}$
其逆函数为：
$P_x^{-1}(u) = (1 - u)^{\frac{1}{1 - \gamma}}$
要从幂律分布 $p_x(x) = (\gamma - 1)x^{-\gamma}$ 中采样 $\xi$，在 $[0, 1]$ 上均匀采样 $u$ 并计算 $\xi = (1 - u)^{\frac{1}{1 - \gamma}}$。

3.4 逆函数法的局限性

累积分布函数 $P_x(x)$ 是 $x$ 的非减函数，其逆 $P_x^{-1}(y)$ 总是存在。但问题是 $p_x(x)$ 可能没有原函数，此时使用逆函数法会更复杂，因为需要求解非线性积分方程来获取每个样本。例如标准高斯分布：
$p_x(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}$
其原函数只能用积分表示，逆函数无法用封闭形式表达。虽然有专门的高斯分布采样算法，但不能用于其他连续概率分布，此时接受 - 拒绝法是一个合适的替代方法。

3.5 流程图：逆函数法采样流程

graph TD;
    A[确定概率密度函数 px(x)] --> B[计算累积分布函数 Px(x)];
    B --> C[求 Px(x) 的逆函数 Px-1(u)];
    C --> D[在 [0, 1] 上均匀采样 u];
    D --> E[计算 ξ = Px-1(u)];
    E --> F[得到按 px(x) 分布的样本 ξ];

4. 接受 - 拒绝法

4.1 原理

当计算累积分布函数 $P_x(x)$ 的逆函数代价高昂，或者 $p_x(x)$ 的原函数无法解析求逆时，接受 - 拒绝法是一种有效的替代方法。

采样具有分布 $p_x(x)$ 的随机变量有明确的几何解释。考虑平面上 $x$ 轴与函数 $p_x(x)$ 之间的区域 $\Omega_p$，其面积为 $1$。在 $\Omega_p$ 中均匀采样一个点，该点横坐标 $x$ 按概率密度 $p_x(x)$ 分布。

为了从 $p_x(x)$ 中采样，我们寻找一个函数 $w(x) \geq p_x(x) \forall x \in R$，其原函数 $W(x)$ 易于求逆。在区域 $\Omega_w$ 中均匀采样一个点 $(x, y)$，若该点在 $\Omega_p$ 内则接受样本，否则拒绝并重复采样，直到接受一个样本。

4.2 实现步骤

具体实现步骤如下：
1. 计算 $A = \int_{-\infty}^{\infty} w(x) dx$，使用逆函数法从 $w(x)/A$ 中采样随机变量 $x$。
2. 在 $[0, w(x)]$ 上均匀采样 $y$。
3. 若 $y < p_x(x)$，则接受采样值 $x$，否则拒绝。

以下是接受 - 拒绝采样的伪代码：

Algorithm 21 acceptance_rejection_sampling()
1: ξ1 ← RAND(0, 1)
2: x = W−1 (Aξ1)
3: ξ2 ← RAND(0, 1)
4: y ← w(x)ξ2
5: if y < px(x) then
6:     return x { accept the sample x and return it}
7: else
8:     goto 1 { reject the sample and repeat }
9: end if

4.3 效率分析

接受 - 拒绝法的效率取决于求逆函数 $W(x)$ 的难易程度以及接受样本的频率。接受样本的比例取决于 $p_x(x)$ 和 $w(x)$ 下方的面积，等于 $1/A$。因此，高效采样要求 $w(x)$ 的原函数 $W(x)$ 易于求逆，且曲线 $w(x)$ 下方的面积 $A$ 尽可能接近 $1$。

4.4 流程图：接受 - 拒绝法采样流程

graph TD;
    A[确定 px(x) 和 w(x)] --> B[计算 A = ∫w(x)dx];
    B --> C[在 [0, 1] 上均匀采样 ξ1];
    C --> D[计算 x = W-1(Aξ1)];
    D --> E[在 [0, 1] 上均匀采样 ξ2];
    E --> F[计算 y = w(x)ξ2];
    F --> G{y < px(x)?};
    G -- 是 --> H[接受样本 x];
    G -- 否 --> C;

5. 离散分布采样

5.1 简单离散采样

在网络科学中，常常需要根据一定的概率质量函数 $p_k$ 从一组整数中采样。假设 $k \in {1, \ldots, M}$ 是一个整数随机数，概率为 ${p_k}$，要根据 ${p_k}$ 采样一个整数 $\kappa$，可以按以下步骤进行：
1. 将区间 $[0, 1]$ 划分为 $M$ 个区间，第 $i$ 个区间的长度等于 $p_i$。
2. 在 $[0, 1]$ 上均匀采样一个随机数 $\xi$。
3. 确定 $\xi$ 所属的子区间 $k$，并将 $k$ 作为采样结果返回。

以下是简单离散采样的伪代码：

Algorithm 22 discrete_sampling()
1: Compute Pi = ∑i m=1 pm, i = 1, ..., M, P0 = 0;
2: find the integer k such that Pk−1 ≤ ξ < Pk.

该采样方案的效率取决于确定样本 $\xi$ 所属子区间的难易程度。在某些特定情况下，如几何分布的随机变量，可以通过解析方法确定该区间。

5.2 几何分布变量采样

对于几何分布的随机变量，概率质量函数为 $p_k = (1 - \tau)\tau^{k - 1}, k = 1, \ldots, \infty$，累积分布函数为 $P_k = \sum_{j = 1}^{k} p_j = 1 - \tau^k$。满足条件 $P_{k - 1} \leq \xi < P_k$ 的整数 $k$ 可以通过解析方法确定为：
$k = \lfloor\frac{\ln(1 - \xi)}{\ln \tau}\rfloor + 1$
因此，算法 22 的第 2 行可以简化为上述公式，并且该公式也适用于 $M = \infty$ 的情况。

5.3 任意概率质量函数采样

对于任意一组概率 ${p_k}$，可以使用二分查找来实现算法 22 的第 2 行检查。具体做法是，计算与 ${p_k}$ 相关的累积概率分布 $P_k$ 的值，并将它们存储在数组 cumul[] 中，其 $M$ 个分量满足 cumul[k] = Pk 。当需要检查样本变量 $\xi$ 对应哪个 $k$ 值时，只需对数组 cumul 使用二分查找算法。这样，算法 22 的第 2 行检查可以在 $O(\log_2 M)$ 操作内完成。

算法 22 的第 1 行复杂度为 $O(M)$，因为需要计算 $M$ 个 $P_i$ 值来确定区间 $k$。通常，需要从给定的 ${p_k}$ 中采样 $N$ 个 $k$ 值，因此算法的第 2 行需要重复 $N$ 次。该过程的总体复杂度为 $O(M) + O(N \log_2 M)$。如果 $M$ 非常大，即 $M \gtrapprox N \log_2 M$，则大部分时间将花费在计算 ${P_k}$ 上。在这种情况下，可以使用更高效的接受 - 拒绝技术。

5.4 离散接受 - 拒绝采样

假设可以有效地计算最大 $p_k$ 的估计值 $\tilde{p} {max}$（$\tilde{p} {max} \geq p_k, k = 1, \ldots, M$），离散接受 - 拒绝技术的步骤如下：
1. 在 $[1, M]$ 上均匀采样一个整数 $k$。
2. 以概率 $p_k / \tilde{p}_{max}$ 接受该样本。

几何上，这相当于在区间 $[0, M\tilde{p}_{max}]$ 上均匀采样一个点。如果该点落在灰色区间内（如图所示），则采样对应的 $k$ 值。

以下是离散接受 - 拒绝采样的伪代码：

Algorithm 23 discrete_sampling_acceptance_rejection()
1: select an integer random number uniformly in [1, ..., M], i.e.:
2:     ξ1 ← RAND(0, 1)
3:     k ← ⌊Mξ1⌋ + 1
4: ξ2 ← RAND(0, 1)
5: if pkξ2 < p then
6:     accept the sample
7: else
8:     reject it and go to 1
9: end if

该过程可以推广到估计值 $p$ 依赖于 $k$ 的情况，即当 $p_k \geq p_k, k = 1, \ldots, M$ 时，这是连续情况下接受 - 拒绝采样的离散版本。此时，需要从 ${p_k / P}$ 中采样，其中 $P = \sum_{k} p_k$，然后如果 $\xi p_k \leq p_k$ 则接受样本。

5.5 流程图：离散接受 - 拒绝采样流程

graph TD;
    A[确定 pk 和 pmax] --> B[在 [0, 1] 上均匀采样 ξ1];
    B --> C[计算 k = ⌊Mξ1⌋ + 1];
    C --> D[在 [0, 1] 上均匀采样 ξ2];
    D --> E{pkξ2 < p?};
    E -- 是 --> F[接受样本 k];
    E -- 否 --> B;

6. 离散幂律分布采样

作为离散分布采样的一个例子，考虑从幂律分布 $p_k = \frac{k^{-\gamma}}{\zeta(\gamma)}, k = 1, \ldots$ 中采样的问题。由于算法 22 需要计算 $\sum_{m = 1}^{i} k^{-\gamma}$，没有简单的表达式，因此使用接受 - 拒绝方法。

具体做法是，考虑图中黑色实线下方的区域，它由两部分组成：$x < 1$ 时的单位面积正方形和 $x > 1$ 时曲线 $y = x^{-\gamma}$ 下方面积为 $A = \frac{1}{\gamma - 1}$ 的区域。为了从灰色区域采样，使用接受 - 拒绝方法：
1. 在黑色实线下方的区域均匀采样一个点 $(x, y)$。
2. 如果采样点在灰色区域内，则接受该样本，否则拒绝并重复该过程。

实际操作中，首先检查点是否在第一个区域（单位面积正方形）内，其概率为 $\frac{1}{1 + A} = \frac{\gamma - 1}{\gamma}$。否则，在第二个区域采样：$x$ 通过逆映射采样，$y$ 在 $[0, x^{-\gamma}]$ 上均匀采样。

以下是离散幂律分布采样的伪代码：

Algorithm 24 power_law_sample()
Input: γ
Output: a sample from the discrete power-law k−γ / ζ(γ)
1: ξ0 ← RAND(0, 1)
2: if ξ0 < (γ - 1) / γ then
3:     k = 1
4: else
5:     ξ1 ← RAND(0, 1)
6:     ξ2 ← RAND(0, 1)
7:     x ← (1 - ξ1)^{\frac{1}{1 - \gamma}}
8:     ktrial ← ⌊x⌋ + 1
9:     y ← ξ2x^{-\gamma}
10:    if y < k_{trial}^{-\gamma} then
11:        return ktrial
12:    else
13:        go to 1
14:    end if
15: end if

6.1 流程图：离散幂律分布采样流程

graph TD;
    A[输入 γ] --> B[在 [0, 1] 上均匀采样 ξ0];
    B --> C{ξ0 < (γ - 1) / γ?};
    C -- 是 --> D[k = 1];
    C -- 否 --> E[在 [0, 1] 上均匀采样 ξ1];
    E --> F[计算 x = (1 - ξ1)^{\frac{1}{1 - \gamma}}];
    F --> G[计算 ktrial = ⌊x⌋ + 1];
    G --> H[在 [0, 1] 上均匀采样 ξ2];
    H --> I[计算 y = ξ2x^{-\gamma}];
    I --> J{y < ktrial^{-\gamma}?};
    J -- 是 --> K[返回 ktrial];
    J -- 否 --> B;

综上所述，随机采样在许多算法中都有重要应用，不同的采样方法适用于不同的概率分布和场景。通过合理选择采样方法，可以提高采样效率和准确性。在实际应用中，需要根据具体问题的特点来选择合适的采样方法，并注意参数的选择对采样结果的影响。