松鼠的窝
首发于松鼠的窝
10809 一种错误的洗牌算法,以及乱排常数 (1)

10809 一种错误的洗牌算法,以及乱排常数 (1)

  「洗牌」,或者说随机打乱一个数组中元素的顺序,是编程中的一个常见需求。标准的洗牌算法是 Fisher-Yates shuffle,用 JavaScript 实现如下:

function shuffle(A) {
    for (var i = A.length - 1; i > 0; i--) {
        var j = Math.floor(Math.random() * (i + 1));
        var t = A[i]; A[i] = A[j]; A[j] = t;
    }
}

其基本思路是,每次从未打乱的部分等可能地选一个元素,把它与未打乱部分的最后一个元素交换。

  Fisher-Yates 洗牌算法的实现十分简单,并且它可以保证均匀性,即元素的各种排列顺序出现的概率都相等。但是,很多人闭门造车地发明了一些「错误」的洗牌算法实现,它们不能保证均匀。例如,最常见的一种错误实现如下:

function shuffle(A) {
   for (var i = 0; i < A.length; i++) {
       var j = Math.floor(Math.random() * A.length);
       var t = A[i]; A[i] = A[j]; A[j] = t;
   }
}

其原理是,在第 i 次循环中,从所有元素中等可能地选一个元素,与第 i 个元素交换。这种算法的错误可以如下证明:对于一个长度为  n 的数组,算法创造了 n^n 个等可能的基本事件,这些事件对应于 n! 种排列顺序。在非平凡情况下, n^n 不能被 n! 整除,所以各种排列顺序不可能等概率。

  另一种错误的洗牌算法是 @Lucas HC这篇答案中指出的:

A.sort(function() {
    return .5 - Math.random();
}); 

JavaScript 中数组自带的 sort 方法允许提供一个比较器,其返回值的正负号代表两个元素的大小关系。在上面的代码中,比较器返回的是 -0.5 到 0.5 之间的一个随机数,也就是说每次比较的结果是随机且均匀的。但是,基于随机比较的整个洗牌算法是不均匀的:它的各种运行结果的概率都形如 2^{-m}m 为算法执行过程中的比较次数),而我们希望每种顺序的概率都是 1/n! ,在非平凡情况下,后者不能由前者通过加法组合出来。@Lucas HC 指出,当 sort 函数采用插入排序的实现时,各个元素都有较大的概率留在初始位置,并通过统计多次运行的结果进行了验证。

  如果你不熟悉编程,我在此可以用大白话把 @Lucas HC 答案中错误算法的流程叙述一下:设想有 n 个人依次来到一个队伍。每个人到来之后,都想向前插队,但前面每一个人放他过去的概率都是 1/2。最终队伍的状态就是洗牌结果。

5 号元素插到了 3 号位置。这个事件要发生,需要队尾的两个元素放它过去,而前面的一个元素不放它过去,故概率为 1/8。

  我对这种算法的结果分布感到好奇,于是计算了一下洗牌后各个元素落在各个位置的概率。用 p^{(n)}_{i,j} 表示总共有 n 个元素时,第 i 个元素洗牌后落在第 j 个位置的概率(这里 i,j 的范围为 1 到 n)。给定 n,各个 p^{(n)}_{i,j} 可以排成一个 n 阶矩阵,记作 P^{(n)} 。这些矩阵可以用递推法计算:

  • n=1 时,显然有 p^{(1)}_{1,1} = 1
  • n>1 时:
    • 首先考虑 n 号元素。它有 1/2 的概率停在 n 号位置,有 1/4 的概率向前插一位停在 n-1 号位置,有 1/8 的概率向前插两位停在 n-2 号位置 …… 一般地,若 n 号元素要停在 j 号位置(2 \le j \le n),它就需要越过 n-j 个元素,并被 j-1 号位置的元素挡住,故概率为 P^{(n)}_{n,j} = 2^{-(n-j+1)},但停在 1 号位置的概率跟停在 2 号位置的概率相同,都是 P^{(n)}_{n,1} = P^{(n)}_{n,2} = 2^{-(n-1)}
    • 然后考虑 i1 \le i < n)号元素。它落在 j 号位置有两种可能:
      • 在前 n-1 个元素洗牌完毕后,i 号元素就已经落在 j 号位置,并且 n 号元素没能越过它。这种情况的概率为 P^{(n-1)}_{i,j} \cdot [1 - 2^{-(n-j)}]
      • 在前 n-1 个元素洗牌完毕后,i 号元素落在 j-1 号位置,但 n 号元素越过了它,使它后移到了 j 号位置。这种情况的概率为 P^{(n-1)}_{i,j-1} \cdot 2^{-(n-j+1)}
    • 于是有递推式:P^{(n)}_{i,j} = P^{(n-1)}_{i,j} \cdot [1 - 2^{-(n-j)}] + P^{(n-1)}_{i,j-1} \cdot 2^{-(n-j+1)}

  上面的递推式用 Matlab 可以如下计算。计算过程使用了大量的矩阵运算,看不懂不必强求。

N = 100;
P = cell(1, N);
P{1} = 1;
for n = 2:N
    x = 0.5 .^ [n-1:-1:1];
    P{n} = zeros(n, n);
    P{n}(1:n-1, 1:n-1) = bsxfun(@times, P{n-1}, 1 - x);
    P{n}(1:n-1, 2:n) = P{n}(1:n-1, 2:n) + bsxfun(@times, P{n-1}, x);
    P{n}(n, 1) = x(1);
    P{n}(n, 2:n) = x;
end

  上述代码算出的 P^{(10)} 如下。左图直接显示了整个矩阵,颜色深浅表示概率大小;右图则是画出了矩阵的每一行,即每个元素最终位置的分布。

使用错误的算法洗牌后,10 个元素位于各个位置的概率

@Lucas HC 的结论,就是说 P^{(n)} 的对角线元素会显著大于 1/n。上图说明, P^{(10)} 的对角线元素都大于 0.2,印证了这个结论。从图中还能看出另外几个事实:

  • P^{(n)} 是对称的,即 i 号元素落在 j 号位置的概率,等于 j 号元素落在 i 号位置的概率;
  • P^{(n)} 的前两行、前两列分别相等,也就是说 1、2 号元素的地位相同,1、2 号位置的地位也相同。

  把 n 增大,可以发现一些更有趣的事情。例如,下图画出的是 P^{(30)}

使用错误的算法洗牌后,30 个元素位于各个位置的概率

从左图可以看到,矩阵中元素沿垂直于对角线的方向迅速衰减,所以「洗牌」后各个元素都倾向于留在初始位置或其附近。而右图则说明,除了开头和结尾几个元素,每个元素最终位置分布曲线的形状都几乎是一样的,它们的最大值(即留在原位的概率)都比 0.2 大一点,且这个值似乎与 n 无关!

  在 Matlab 中执行 format long 命令,可以让计算结果具有 15 位有效数字。用此方法可以观察到 P^{(n)} 的对角线中部元素趋向于一个常数 C ,它约等于 0.220643036096533。尽管 P^{(n)} 的绝大多数对角元都只比常数 C 大一丁点儿,不过如果要在矮子里面再挑矮子,那么最小的对角元大约出现在第 \sqrt{n} 个位置。当 n=64 时,15 位有效数字就已经无法区分 P^{(64)}_{8,8} 与常数 C 了。

  网上并没有找到对常数 C \approx 0.220643036096533 的讨论,所以我就暂且给它起个名字,叫「乱排常数」。这个常数是有理数还是无理数呢?它有什么样的数学表达式?这些有趣的问题,就留到下一篇中讨论。

编辑于 2017-12-02

文章被以下专栏收录