抽样问题——《编程珠玑》读书笔记

问题：输入两个整数m和n，并且m<n。输出一个由m个随机数字组成的有序列表，这些随机数的范围是[0, n-1]，并且每个整数最多出现一次。

方法一：

Knuth著作《Seminumerical Algorithms》中提出的方法，顺序遍历n个数，通过随机测试条件的元素被选择。

以一个例子来解释所说的随机测试条件，比如m=2，n=5。第一个元素0被选择的概率是2/5；第二个元素1被选择的概率取决于第一个元素有没有被选择，如果0被选择，则1被选择的概率为1/4，否则为2/4，所有1被选择的概率为（2/5）*（1/4）+（3/5）*（2/4）=2/5；同理第三个元素2被选择的概率取决于前两个的选择情况，如果都没被选择，则2被选择的概率为2/3，如果前两个有一个被选择，则2被选择的概率为1/3，如果前两个都被选择，则2被选择的概率为0，故2被选择的概率为（3/5）*（3/5）*（2/3）+2*（2/5）*（3/5）*（1/3）=2/5。依次类推，每个元素被选择的概率都为2/5。

总的来说，从剩下的r个元素中选择s个元素，那么下一个元素被选中的概率为s/r，从整个数据集合角度来讲，每个元素被选择的概率都是相同的。

这个思想的为代码如下：

select = m
remaining = n
for i = [0, n)
        if (bigrand() % remaining) < select
                print i
                --select
        --remaining

首先，该算法可以保证有m个元素被选中，不会多也不会少。证明如下，首先证明不会多于m个：因为select等于0时不能选择更多的整数；再证明不会少于m个：当select/remaining=1时，总会选中一个元素，因为bigrand()%remaining<remaining总成立，所以i总会被选中。

其次，每个元素被选择的概率是相等的，均为m/n，证明如上举例证明所示。 C++实现代码如下，同时算出从268435455（约2.7亿，int可以表示的最大整数除以8，本来准备拿int的最大整数约21.5亿测试，但方法三要先new出这么大的空间，超出程序可以分配的最大堆栈空间）个数中选出10万个整数，测试该方法所用的时间，便于与后面的方法进行性能比较。

#include <iostream>
#include <ctime>
#include <cstdlib>
#include <limits>

using namespace std;

void genknuth(int m, int n)
{
	time_t t_start, t_end;
	t_start = time(NULL);
	
	for (int i = 0; i != n; ++i)
		if ((rand() % (n-i)) < m)
		{
			cout << i << " ";
			--m;
		}

	cout << endl;
	t_end = time(NULL);
	cout << "collapse time: " << difftime(t_end, t_start) << " s" << endl;
}

int main()
{
	int m = 100000;
	int n = numeric_limits<int>::max() / 8;
	srand(time(NULL));
	genknuth(m, n);
	cout << "n = " << n << endl;
	return 0;
}

该算法的空间复杂度为O（m），时间复杂度为O（n），用这算法从2.7亿个数中随机找出10万个数所用时间为4秒。

方法二：

方法一所需时间和搜索空间成正比，有些应用仍不能接受，因此需要继续改进。其中一种方法是随机插数据到一个容量为m的集合中。为代码如下所示：

initialize set S to empty
size = 0
while size < m do
        t = bigrand() % n
        if t is not in S
                insert t into S
                ++size
print the elements of S in sorted order

C++代码实现如下所示，集合S的实现采用stl提供的set，底层用红黑树实现，不可重复插入相同的数据，当要插入的数据在set中已经存在时，则插入无效，数据不会被插入集合中，插入的时间复杂度为O（logm）：

#include <iostream>
#include <ctime>
#include <cstdlib>
#include <limits>
#include <set>

using namespace std;

void gensets(int m, int n)
{
	time_t t_start, t_end;
	t_start = time(NULL);

	set<int> S;
	while (S.size() < m)
		S.insert(rand() % n);
	for (set<int>::iterator iter = S.begin(); iter != S.end(); ++iter)
		cout << *iter << " ";

	cout << endl;
	t_end = time(NULL);
	cout << "collapse time: " << difftime(t_end, t_start) << " s" << endl;
}

int main()
{
	int m = 100000;
	int n = numeric_limits<int>::max() / 8;
	srand(time(NULL));
	gensets(m, n);
	return 0;
}

算法的时间复杂度为O（mlogm），空间复杂度为O（m）。同样从2.7亿数据范围内选10万个数，所花的时间为2秒，可见速度会比原来的Knuth方法快。

方法三：

弄乱一个n个元素的数组，然后排序输出前m个元素。后来Ashley
Shepherd和Alex Woronow发现，只需弄乱数组前m个元素，关于产生随机序列的方法可以参考我的文章《洗牌程序》和维基百科。

本方法的C++代码实现如下：

#include <iostream>
#include <ctime>
#include <cstdlib>
#include <limits>
#include <algorithm>

using namespace std;

// generate a random number between i and j,
// both i and j are include.
int randint(int i, int j)
{
	int ret = i + rand() % (j - i + 1);
	return ret;
}

void genshuf(int m, int n)
{
	time_t t_start, t_end;
	t_start = time(NULL);

	int i, j;
	int *x = new int
;
	for (i = 0; i != n; ++i)
		x[i] = i;
	for (i = 0; i != m; ++i)
	{
		j = randint(i, n-1);
		int t = x[i]; x[i] = x[j]; x[j] = t; // swap x[i] and x[j]
	}
	sort(x, x + m);
	for (i = 0; i != m; ++i)
		cout << x[i] << " ";
	
	cout << endl;
	t_end = time(NULL);
	cout << "collapse time: " << difftime(t_end, t_start) << " s" << endl;

	delete []x;
	x = NULL;
}

int main()
{
	int m = 100000;
	int n = numeric_limits<int>::max() / 8;
	srand(time(NULL));
	genshuf(m, n);
	return 0;
}

算法的时间复杂度是O（n+mlogm），空间复杂度是O（n）。同样从数据范围内选10万个数，所花时间为4秒，时间和方法一差不多，其中有一份部分时间花在初始化数组上，如果采用《编程珠玑》问题1.9的方法，当用到某个数时才初始化，这样算法的时间复杂度可以减少到O（mlogm）。不过空间复杂度O（n）还是太大。

关于具体采用方法二还是方法三，stackflow上有个大牛用数学的方法证明了一下，当m<<n时，采用方法二会比三性能好。

参考文章：

http://www.cnblogs.com/2010Freeze/archive/2012/02/27/2370284.html

http://hi.baidu.com/23star/blog/item/47f7314e5c3b0e01b2de0574.html

Taking Random Samples

A Sample of Brilliance,
Programming Perls