随机生成排列组合

假设 random(a, b) 以相同概率随机返回闭区间 [a, b] 内的整数.

均匀随机排列 (洗牌算法), O(n) 如果 1~n 的每一种排列都以 1/n! 的概率出现, 则称该算法实现了均匀随机排列. 注意, 只证明每个元素被选中概率为 1/n 是不充分的. 算法的思路是为每个元素随机生成唯一关键字, 再按关键字排序. 时间复杂度取决于排序算法. 一般的排序算法为 O(n log n), 也有特殊算法可以在 O(n) 时间内对 1 ~ n3 的整数排序. 

shuffle_by_sort(L):
    n = L.len
    sort(L, key=[random(1, n**3) for i = 0..n-1])

设 `n` 个关键字为 `k_1, cdots, k_n`. 用 `A_i` 表示 L[i+1] 分配到 `k_i` 这一事件. 则 `P(nnn_(i=1)^n A_i)` `= P(A_1) P(A_2 | A_1) cdots P(A_n | nnn_(i=1)^(n-1) A_i)` `= 1/n 1/(n-1) cdots 1/2 1/1 = 1/(n!)`. 即每一特定排列出现的概率都是 `1//n!`.

均匀随机 k 排列, O(k) 从 1~n 选出 k 个元素, 它们的任一排列称为一个 k 排列. 该算法在线性时间 O(k) 内实现了均匀随机 k 排列, 每个可能的 k 排列出现的概率都是 (n-k)! / n!. 其中每一步迭代时将下标为 i 的元素与下标 ≥ i 的元素随机交换. 最后生成的 k 排列保存在数组的前 k 个位置中. 特别取 k = n (或 k = n-1), 就得到均匀随机排列. 

shuffle(L, k):
    n = L.len
    for i = 0..k-1:
        swap(L, i, random(i, n-1))

只需证: 当第 `i` 次循环结束时, 子数组 [0..i] 包含一个均匀随机 `i+1` 排列.
`i = 0` 时, 我们将 L[0] 与数组的任意元素随机交换, 结论成立. 设结论对 `i-1` 成立, 来证 `i` 的情形. 任取一个 `i+1` 排列 `x_0, x_1, cdots, x_i`, 记子数组 [0..i] 恰好包含这个排列的事件为 `A_i`. 子数组 [0..i-1] 恰好包含排列 `x_0, x_1, cdots, x_(i-1)` 的事件为 `A_(i-1)`. 由条件概率公式, `P(A_i) = P(A_(i-1)) P(A_i|A_(i-1))` `= ((n-i)!)/(n!) 1/(n-i)` `= ((n-i-1)!)/(n!)`. 最后, 令 `i = k-1` 即得算法的正确性.

Fisher-Yates 洗牌算法 准备一个空队列, 从 `n` 个数字中随机选一个加到队尾. 再从剩下的 `n-1` 个数字中随机选一个加到队尾, 依此类推. 

fisher_yates(arr):
    i = arr.length
    if i < 2:
        return
    while --i:
        j = rand() % (i+1)
        swap(arr, i, j)

希望借助下例, 扫除一些盲区:

    从 1~n 中随机选取 k 个数 (无放回),
  1. 若每个 k 子集被选中的概率相等, 则每个数被选中的概率相等, 反之不对;
  2. 若每个数被选中的概率都相等, 则这个概率等于 k/n.
  1. 设每个 `k` 子集被选中的概率相等, 由于对任意数字 `n_i`, 含有 `n_i` 的子集数为 `(n-1; k-1)`, 从而 `n_i` 被选中的概率为 `(n-1;k-1) // (n;k) = k/n`. 反之, 从 `1~4` 中以相等概率 `1//4` 选取子集 `{1, 2}`, `{2, 3}`, `{3, 4}`, `{4, 1}`, 则每个数被选中的概率都等于 `1//2`, 但 `{1, 3}`, `{2, 4}` 的组合并未出现, 所以并不是每个子集被选中的概率都相等.
  2. 设每个数被选中的概率都等于 `x_0`. 记 `m = (n;k)`, `m` 个 `k` 子集及其出现的概率分别为 `S_j` 和 `x_j`, `j = 1, cdots, m`, 则 `{ sum_(i=1)^m x_i = 1; sum_(i in S_j) x_i = x_0", "j = 1, cdots, m; :}` 后 `m` 个方程相加得 `(n-1;k-1) sum_(i=1)^m x_i = x_0 (n;k)`. 即得 `x_0 = k//n`.

Knuth 抽样算法, O(n) 从 1~n 中随机抽取 k 个数字, 使得每个数被选中的概率都等于 k/n. 上文的均匀随机 k 排列算法可以解决此问题. 现假设数据只能顺序读取, 不能通过下标随机访问, 如何实现随机抽样?
我们取红球 k 个, 黑球 n-k 个, 做 n 次不放回摸球试验. 若第 i 次摸出红球, 则选中数字 i, 否则不选中它. 由于抽签概率与顺序无关, 每个数字被选中的概率都相等. 每个 k 子集出现的概率相等吗? 

sample(L, k):
    n = L.len
    for i = 0 to n-1:
        if random(1, n-i) ≤ k:
            print(L[i])
            --k

蓄水池抽样算法, O(n) 现假设数据是流式的, 只能顺序读取, 而且总数 n 未知. 为实现随机抽样 (使得每个数被选中的概率都是 k/n), 我们维护一个大小为 k 的蓄水池, 初始时装着前 k 个元素. 从 i = k 开始, 每次循环随机生成 [0, i] 中的下标 d, 若 d 落在蓄水池范围内, 则该元素将取代蓄水池中的对应元素. 

reservoir(L, k):
    buf = [L[i] for i = 0..k-1]
    for (i = k; L[i] != null; ++i):
        d = random(0, i)
        if d < k:
            buf[d] = L[i]

注意每个元素只进入蓄水池一次, 而且只有后面的元素能替换前面的元素, 我们得到: 每个元素被选取的概率等于它晋级的概率乘以它不被替换的概率. 于是 L[i] 被选取的概率等于 `{ 1 * k/(k+1) (k+1)/(k+2) cdots (n-k)/n = k/n, if i lt k; k/i * i/(i+1) (i+1)/(i+2) cdots (n-k)/n = k/n, if i ge k :}` 每个 k 子集出现的概率相等吗?

并行蓄水池算法 [简书 邱simple]
令 `t` 个线程分别用蓄水池算法处理 `n_i` 份数据, 每个线程各得到 `k` 份样本, 其中 `sum_(i=1)^t n_i = n`, `n_i ge k`. 作累加 `m_i = sum_(j=1)^i n_j`. 随机选取 `d in [0, n)`, 若 `d in [m_(i-1), m_i)`, 则在第 `i` 个线程的蓄水池中等概率不放回地选取一个数据, 如此重复 `k` 次, 最终得到 `k` 个样本.

第 `i` 个线程中, 数据被选取的概率为 `k//n_i`. 从这个线程的蓄水池中选取数据, 放进最终蓄水池的概率为 `n_i//n`. 相乘即得到每个数据被选中的概率 `k/n_i * n_i/n * 1/k * k = k/n`.

伪随机数

线性同余法生成伪随机数 (D. H. 莱默, 1949)
选取整数 `2 le a lt m` 和 `0 le c lt m`, 种子 `0 le x_0 lt m`, 则递推公式 `x_(n+1) = a x_n + c (mod m)` 确定了一组伪随机数 `{x_n}`. 注意一旦有 `x_i = x_j`, 就有 `x_(i+1) = x_(j+1)`, 序列将会重复. 我们把序列出现重复前的最大长度称为它的周期. 由鸽巢原理, 周期不超过 `m`. 给出周期恰好为 `m` 的充要条件. 特别, 当 `c = 0` 时, 称为纯乘性同余方法. 下面是一个实用的递推公式 `x_(n+1) = 252,246,292 x_n (mod 2^31-1)`. 密码学中要求伪随机数不可预测. 线性同余方法不能用于密码学, 因为已知连续的若干项就可以求得其它项.

    线性同余生成器产生的伪随机序列周期为 `m` 当且仅当以下三条同时满足:
  1. `(c, m) = 1`;
  2. 若素数 `p | m`, 则 `a -= 1 (mod p)`;
  3. 若 `4 | m`, 则 `a -= 1 (mod 4)`.

正态随机数

Box-Muller 算法, 每次生成一对正态随机数 设 `X, Y overset"iid"~ N(0, 1)`, 其联合密度为 `p(x, y) = 1/(2pi) "e"^(-(x^2+y^2)/2)`. 作变量替换 `X = r cos theta`, `Y = r sin theta`, 则 `1 = int_(-oo)^oo int_(-oo)^oo 1/(2pi) "e"^(-(x^2+y^2)/2) dx dy` `= int_0^(2pi) int_0^oo 1/(2pi) "e"^(-r^2/2) r "d"r "d"theta`. 从而得到 `r, theta` 的分布函数 `P(r le r_0) = int_0^(2pi) int_0^(r_0) 1/(2pi) "e"^(-r^2/2) r "d"r "d"theta` `= 1 - "e"^(-r_0^2/2)` `quad (r_0 ge 0),
`P(theta le theta_0) = int_0^(theta_0) int_0^oo 1/(2pi) "e"^(-r^2/2) r "d"r "d"theta` `= theta_0/(2pi)` `quad (0 le theta_0 le 2pi)`. 如何生成满足上面这两个分布的随机变量 `r, theta`? 注意到 `theta ~ U(0, 2pi)`; 至于 `r`, 可以利用反函数: 令 `z = 1 - "e"^(-r^2/2)`, 反解得 `r = sqrt(-2 ln(1-z))`. 于是当 `z ~ U(0, 1)` 时, `r` 的分布函数为 `1 - "e"^(-r^2/2)`. 总之, 当 `U_1, U_2 ~ U(0, 1)` 时, 令 `X = r cos theta = sqrt(-2 ln U_1) cos 2pi U_2`,
`Y = r sin theta = sqrt(-2 ln U_1) sin 2pi U_2`, 则 `x, y` 服从标准正态分布.

下面的代码使用 Box-Muller 算法生成一对正态随机数; 该函数直接返回其中一个随机数, 另一个被缓存, 等下次调用时返回.