引言
在统计学中,简单随机取样(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# 程序如下所示:
1 public static IEnumerable<int> Sample01(int m, int n) 2 { 3 var s = new HashSet<int>(); 4 var r = new Random(); 5 while (s.Count < m) s.Add(r.Next(n) + 1); 6 return s; 7 }
但是算法S有一个缺陷:当 M = N 时,集合 S 缺一个整数。算法闭着眼睛乱猜整数,直到偶然碰上正确的那个为止,这平均需要猜 N 个随机数。这个分析假设随机数发生器是真正随机的。对于某些非随机序列,这个算法甚至不会停止。
算法S2
算法S的缺陷是可以避免的,只要发现 M > N / 2,就令 M := N - M,执行算法S,然后再从 1 .. N 中扣除集合 S 就行了,如下所示:
1 public static IEnumerable<int> Sample02(int m, int n) 2 { 3 if (m <= n / 2) return Sample01(m, n); 4 return Enumerable.Range(1, n).Except(Sample01(n - m, n)); 5 }
算法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# 程序如下所示:
1 public static IEnumerable<int> Sample03(int m, int n) 2 { 3 return Sample03A(new Random(), m, n); 4 } 5 6 static HashSet<int> Sample03A(Random r, int m, int n) 7 { 8 if (m == 0) return new HashSet<int>(); 9 var s = Sample03A(r, m - 1, n- 1); 10 if (!s.Add(r.Next(n) + 1)) s.Add(n); 11 return s; 12 }
注意 ,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# 程序如下所示:
1 public static IEnumerable<int> Sample04(int m, int n) 2 { 3 var s = new HashSet<int>(); 4 var r = new Random(); 5 for (var j = n - m + 1; j <= n; j++) 6 if (!s.Add(r.Next(j) + 1)) s.Add(j); 7 return s; 8 }
《编程珠玑(续)》还给出了该算法的 awk 程序(这本书中的程序示例都用是 C 语言和 awk 语言写的):
1 BEGIN { m = ARGV[1]; n = ARGV[2] 2 for (j = n-m+1; j <= n; j++) { 3 t = 1 + int(j * rand()) 4 if (t in s) s[j] = 1 5 else s[t] = 1 6 } 7 for (i in s) print i 8 }
算法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# 程序如下所示:
1 public static IEnumerable<int> Sample05(int m, int n) 2 { 3 var s = new List<int>(m); 4 var r = new Random(); 5 for (var j = n - m + 1; j <= n; j++) 6 { 7 var t = r.Next(j) + 1; 8 var i = s.IndexOf(t); 9 s.Insert(i + 1, (i == -1) ? t : j); 10 } 11 return s; 12 }
运行实例
上述算法在 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的效果是一样的,只不过一个是递归的算法,一个是迭代的算法。
测试程序
产生上述运行实例的测试程序如下所示:
1 using System; 2 using Skyiv.Libs; 3 4 namespace Skyiv.Tester 5 { 6 static class SampleTester 7 { 8 static void Main(string[] args) 9 { 10 var m = (args.Length > 0) ? int.Parse(args[0]) : 10; 11 var n = (args.Length > 1) ? int.Parse(args[1]) : m * 10; 12 foreach (var x in SimpleRandomSample.Sample05(m, n)) Console.Write(x + " "); 13 Console.WriteLine(); 14 } 15 } 16 }
上述程序中第 12 行的 Sample05 需要根据算法的不同代之以不同的方法名称。
1 using System; 2 using System.Linq; 3 using System.Collections.Generic; 4 5 namespace Skyiv.Libs 6 { 7 static class SimpleRandomSample 8 { 9 // ... 10 } 11 }
上述程序中第 9 行处需要根据不同的算法填入相应的方法。
实际应用
在我提供技术支持的综合业务技术竞赛中,有一个项目是综合业务基本知识,是上机比赛的。选手们使用 IE 浏览器登录比赛网站,在线答题。共 250 道题,其中判断题和单项选择题各 100 道,多项选择题 50 道,比赛时间限制为三十分钟。业务部门事先提供判断题和单项选择题各 1000 道,多项选择题 500 道,用于建立题库。各个选手的题目都是使用上述的简单随机取样算法从题库中抽取的,每个选手抽到的题目都不相同。
注:本文中的算法S、算法F1、算法F2、算法P来源于《编程珠玑(续)》第13章 绝妙的取样。