蓄水池抽样问题

蓄水池抽样问题

编程珠玑中有这么一个题目：

     How could you select one of n objects at random, where you see the objects sequentially but you do not know the value of n beforehand? For concreteness, how would you read a text file, and select and print one random line, when you don’t know the number of lines in advance?

问题可以简化如下：在不知道文本总行数的情况下，如何从文件中随机的抽取一行。再简单的说就是在不知道n的大小的情况下，如何从n个数据中随机抽取一个数。

     首先我们想一个在n知道的情况下，如何随机选择一个数。很简单，通过rand函数来随机产生一个(0,n-1)之间的一个数x，然后取第x个数即可。但是，当前的情况是不知道n的大小，如何保证能够随机抽取呢？从另外一个方面考虑，随机抽取一个数，也就是说只要能够保证每个数能够以相等的概率被抽中就可以了。

蓄水池抽样问题只不过是上面情况的一个扩展，简单说就是要求从n个数据中随机抽取k个数，其中n的大小是未知的。这就是蓄水池抽样算法。

蓄水池抽样算法伪代码如下：

source[n];         //原始数组，n个元素，n未知，数组下标从0开始
array[k];          //用来存放抽取结果
interger i,j;

//填充数组array
for each i in 0 to k-1 do
     array[i] = source[i];
done;

//以逐渐减小的概率替换array中的数据元素
for each i in k to n do
     j = rand(0,i);          //随机生成一个（0,i）之间的随机数
     if j < k then
          array[j] = source[i];     //替换array中第j个元素
     fi
done

算法首先创建一个长度为k的数组（蓄水池）用来存放结果，并将其初始化原始数组中前k个元素，然后从k+1个元素开始迭代直到数组结束，在迭代到元素数组的第i个元素，随机生成一个（0，i）之间的随机数j，如果j<k，那么蓄水池第j个元素被替换为元素数组中第i个元素。

蓄水池算法的证明：

定理：该算法保证每个元素都以k/n的概率被选入蓄水池数组中。

证明如下：归纳法

首先对于任意的i，第i个元素进入蓄水池的概率为k/i；而在蓄水池内每个元素被替换掉的概率为1/k；因此在第i轮第j个元素被替换掉的概率为(k/i)*(1/k)=1/i。接下来用数学归纳法证明当循环结束时，每个元素进入蓄水池的概率为k/n。

假设在i-1次迭代后，任意一个元素进入蓄水池的概率为k/(i-1)。根据上面的结论，在第i次迭代时，该元素被替换的概率为1/i；那么其不被替换的概率为1-1/i = (i-1)/i；在第i次迭代后，该元素在蓄水池内的概率为k/(i-1)*(i-1)/i = k/i。证明结束。

下面以简单的例子实现：从10个数中随机一次抽取6个数，抽取10000次，统计每个元素被抽取的次数。

#include <iostream>
#include <vector>
using namespace std;

/*
*	蓄水池抽样算法
*	v：原始数据数组
*	n：总元素个数
*	k：一次抽取的元素个数
*	return：一次抽取的结果
*/
vector<int> ReservoirSampling(vector<int> v, int n, int k)
{
	vector<int> reservoirArray(v.begin(), v.begin() + k);		//初始化蓄水池数组为元素数组的前k个元素
	int len = v.size();
	if (len != n || k >= n)
		return reservoirArray;
	int i = 0;
	int j = 0;
	for (i = k; i < n; ++i) {				//从第k+1个元素开始迭代
		j = rand() % (i+1) ;				//随机产生一个(0,i)之间的随机数
		if (j < k)
			reservoirArray[j] = v[i];		//替换
	}
	return reservoirArray;
}

int main()
{
	vector<int> v = { 0,1,2,3,4,5,6,7,8,9 };
	const int iter_count = 10000;			//抽取总次数
	int count[10] = { 0 };						//记录每个元素被抽取的次数
	for (int i = 0; i < iter_count; ++i) {
		vector<int> res = ReservoirSampling(v, 10, 6);		//进行抽取
		for (auto x : res)
			count[x]++;
	}
	//输出每个元素被抽中的次数，大约为k/n*iter_count
	for (int i = 0; i < 10; i++)
		cout << "数字" << i << "被抽取次数为：" << count[i] << endl;
}

运行结果如下：

可以发现：基本上每个元素抽取的次数约为6000次，达到了随机抽取的效果。