绝妙的取样

引言

在统计学中，简单随机取样(simple random sample)，是从总体 N 单位中随机地抽取 M 个单位作为样本，使得每一个样本都有相同的概率被抽中。现在让我们要讨论产生 1 .. N 范围内 M 个不重复的随机整数的算法。

算法S

首先，我们来看一下算法S：

算法S
initialize set S to empty
Size := 0
while Size < M do
    T := RandInt(1, N)
    if T is not in S then
        insert T in S
        Size := Size + 1

这个算法中，RandInt(L, U) 返回在 L .. U 上均匀分布的一个整数。这个算法非常简洁，很容易理解。循环不变式是：S 总是包含 1 .. N 范围内 Size 个整数的随机样本。

相应的 C# 程序如下所示：

    public static IEnumerable<int> Sample01(int m, int n)
    {
      var s = new HashSet<int>();
      var r = new Random();
      while (s.Count < m) s.Add(r.Next(n) + 1);
      return s;
    }

但是算法S有一个缺陷：当 M = N 时，集合 S 缺一个整数。算法闭着眼睛乱猜整数，直到偶然碰上正确的那个为止，这平均需要猜 N 个随机数。这个分析假设随机数发生器是真正随机的。对于某些非随机序列，这个算法甚至不会停止。

算法S2

算法S的缺陷是可以避免的，只要发现 M > N / 2，就令 M := N - M，执行算法S，然后再从 1 .. N 中扣除集合 S 就行了，如下所示：

    public static IEnumerable<int> Sample02(int m, int n)
    {
      if (m <= n / 2) return Sample01(m, n);
      return Enumerable.Range(1, n).Except(Sample01(n - m, n));
    }

算法F1

Bob Floyd 发现了一个绝妙的取样算法，对于 S 中每个随机数只恰好调用一次 RandInt：

算法F1
function Sample(M, N)
    if M = 0 then
        return the empty set
    else
        S := Sample(M-1, N-1)
        T := RandInt(1, N)
        if T is not in S then
            insert T in S
        else
            insert N in S
    return S

这个递归的算法很容易理解：为了从 1 .. N 中产生 M 个元素的样本，首先从 1 .. N-1 中产生 M-1 个元素的样本，然后加上第 M 个元素。

注意，在这个算法中，如果 T 等于 N 的话，T 肯定不在集合 S 中，不会出现因为 T 在集合 S 中而执行 insert N in S 语句导致失败的情况 。

请思考一下，如果把 T := RandInt(1, N) 改为 T := RandInt(1, N-1) 的话，对这个算法有何影响。

相应的 C# 程序如下所示：

    public static IEnumerable<int> Sample03(int m, int n)
    {
      return Sample03A(new Random(), m, n);
    }

    static HashSet<int> Sample03A(Random r, int m, int n)
    {
      if (m == 0) return new HashSet<int>();
      var s = Sample03A(r, m - 1, n- 1);
      if (!s.Add(r.Next(n) + 1)) s.Add(n);
      return s;
    }

注意 ，HashSet.Add 方法在被插入的元素已经存在时并不报错，而是返回 false。

算法F2

可以通过引入一个新变量 J 将递归的算法F1改写为迭代形式，这就是算法F2：

算法F2
initialize set S to empty
for J := N - M + 1 to N do
    T := RandInt(1, J)
    if T is not in S then
        insert T in S
    else
        insert J in S

这个算法比算法S更有效率，也很简洁。相应的 C# 程序如下所示：

    public static IEnumerable<int> Sample04(int m, int n)
    {
      var s = new HashSet<int>();
      var r = new Random();
      for (var j = n - m + 1; j <= n; j++)
        if (!s.Add(r.Next(j) + 1)) s.Add(j);
      return s;
    }

《编程珠玑(续)》还给出了该算法的 awk 程序(这本书中的程序示例都用是 C 语言和 awk 语言写的)：

BEGIN { m = ARGV[1]; n = ARGV[2]
    for (j = n-m+1; j <= n; j++) {
        t = 1 + int(j * rand())
        if (t in s) s[j] = 1
        else s[t] = 1
    }
    for (i in s) print i
}

算法P

前面的算法都不能保证样本的元素以随机的顺序出现，因为集合的元素是无序的，程序可以使用任意顺序输出集合中的元素。下面的算法保证了样本中的元素以随机的顺序出现：

算法P
initialize sequence S to empty
for J := N - M + 1 to N do
    T := RandInt(1, J)
    if T is not in S then
        prefix T to S
    else
        insert J in S after T

这个算法和算法F2类似。为了产生 1 .. N 范围内的一组 M 元排列，它会先从 1 .. N-1 中产生一组 M-1 元排列。但是主要数据结构是序列而非集合。

算法P以等概率生成 1 .. N 范围的每一个 M 元排列。循环不变式是：第 i 轮循环后，J = N-M+i 且 S 可能是 1 .. J 中 i 个不同整数的任意排列，并且只有一种途径可以生成这个排列。

相应的 C# 程序如下所示：

    public static IEnumerable<int> Sample05(int m, int n)
    {
      var s = new List<int>(m);
      var r = new Random();
      for (var j = n - m + 1; j <= n; j++)
      {
        var t = r.Next(j) + 1;
        var i = s.IndexOf(t);
        s.Insert(i + 1, (i == -1) ? t : j);
      }
      return s;
    }

运行实例

上述算法在 M = 10 和 N = 100 时的典型运行结果如下所示：

5 17 52 62 93 75 96 87 73 44   算法S
32 34 27 76 38 72 80 12 25 18  算法S2
9 84 88 24 78 91 47 56 8 46    算法F1
54 75 43 32 15 48 62 45 41 100 算法F2
38 59 17 45 20 5 37 79 36 95   算法P

上述算法在 M = 10 和 N = 19 时的典型运行结果如下所示：

13 7 15 12 11 5 3 14 17 10     算法S
1 2 4 7 8 9 10 11 17 19        算法S2
7 10 3 13 8 11 5 15 17 14      算法F1
8 1 10 3 11 13 16 4 18 2       算法F2
9 17 18 6 12 5 7 13 14 11      算法P

上述算法在 M = 10 和 N = 10 时的典型运行结果如下所示：

9 6 7 10 3 2 8 1 4 5           算法S
1 2 3 4 5 6 7 8 9 10           算法S2
1 2 3 4 5 6 7 8 9 10           算法F1
1 2 3 4 5 6 7 8 9 10           算法F2
3 9 5 6 10 1 2 4 8 7           算法P

实际上，算法F1和算法F2的效果是一样的，只不过一个是递归的算法，一个是迭代的算法。

测试程序

产生上述运行实例的测试程序如下所示：

using System;
using Skyiv.Libs;

namespace Skyiv.Tester
{
  static class SampleTester
  {
    static void Main(string[] args)
    {
      var m = (args.Length > 0) ? int.Parse(args[0]) : 10;
      var n = (args.Length > 1) ? int.Parse(args[1]) : m * 10;
      foreach (var x in SimpleRandomSample.Sample05(m, n)) Console.Write(x + " ");
      Console.WriteLine();
    }
  }
}

上述程序中第 12 行的 Sample05 需要根据算法的不同代之以不同的方法名称。

using System;
using System.Linq;
using System.Collections.Generic;

namespace Skyiv.Libs
{
  static class SimpleRandomSample
  {
    // ...
  }
}

上述程序中第 9 行处需要根据不同的算法填入相应的方法。

实际应用

在我提供技术支持的综合业务技术竞赛中，有一个项目是综合业务基本知识，是上机比赛的。选手们使用 IE 浏览器登录比赛网站，在线答题。共 250 道题，其中判断题和单项选择题各 100 道，多项选择题 50 道，比赛时间限制为三十分钟。业务部门事先提供判断题和单项选择题各 1000 道，多项选择题 500 道，用于建立题库。各个选手的题目都是使用上述的简单随机取样算法从题库中抽取的，每个选手抽到的题目都不相同。

---

注：本文中的算法S、算法F1、算法F2、算法P来源于《编程珠玑(续)》第13章 绝妙的取样。

参考资料

1. 编程珠玑(续)，[美] Jon Bentley 著，钱丽艳 刘田 等译，人民邮电出版社，2011年5月第2版
2. Wikipedia: Simple random sample
3. MSDN: Random 类 (System)
4. MSDN: HashSet(T) 类 (System.Collections.Generic)
5. MSDN: List(T) 类 (System.Collections.Generic)