蓄水池抽样问题

落日映苍穹つ 2022-06-11 07:23 192阅读 0赞

**蓄水池抽样问题**

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

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;
    }

运行结果如下：

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

[Center]: /images/20220611/d5a6f5d4a65740dbab18162c0e8f913d.png