水塘抽样(Reservoir Sampling)问题

水塘抽样是一系列的随机算法，其目的在于从包含n个项目的集合S中选取k个样本，其中n为一很大或未知的数量，尤其适用于不能把所有n个项目都存放到主内存的情况。

在高德纳的计算机程序设计艺术中，有如下问题：可否在一未知大小的集合中，随机取出一元素？。或者是Google面试题： I have a linked list of numbers of length N. N is very large and I don’t know in advance the exact value of N. How can I most efficiently write a function that will return k completely random numbers from the list（中文简化的意思就是：在不知道文件总行数的情况下，如何从文件中随机的抽取一行？）。两题的核心意思都是在总数不知道的情况下如何等概率地从中抽取一行？即是说如果最后发现文字档共有N行，则每一行被抽取的概率均为1/N？

我们可以：定义取出的行号为choice，第一次直接以第一行作为取出行 choice ，而后第二次以二分之一概率决定是否用第二行替换 choice ，第三次以三分之一的概率决定是否以第三行替换 choice ……，以此类推。由上面的分析我们可以得出结论，在取第n个数据的时候，我们生成一个0到1的随机数p，如果p小于1/n，保留第n个数。大于1/n，继续保留前面的数。直到数据流结束，返回此数，算法结束。

 

问题一

首先考虑k为1的情况，即：给定一个长度很大或者长度未知数据流，限定对每个元素只能访问一次，写出一个随机选择算法，使得所有元素被选中的概率相等。

设当前读取的是第n个元素，采用归纳法分析如下：

1. n = 1 时，只有一个元素，直接返回即可，概率为1。
2. n = 2 时，需要等概率返回前两个元素，显然概率为1/2。可以生成一个0～1之间的随机数p，p < 0.5 时返回第一个，否则返回第二个。
3. n = 3 时，要求每个元素返回的概率为1/3。注意此时前两个元素留下来的概率均为1/2。做法是：生成一个0～1之间的随机数，若<1/3，则返回第三个，否则返回上一步留下的那个。元素1和2留下的概率均为：1/2 * (1 - 1/3) = 1/3，即上一步留下的概率乘以这一步留下（即元素3不留下）的概率。
4. 假设 n = m 时，前n个元素留下的概率均为：1/n = 1/m；
5. 那么 n = m+1 时，生成0～1之间的随机数并判断是否<1/(m+1)，若是则留下元素m+1，否则留下上一步留下的元素。这样一来，元素m+1留下的概率为1/(m+1)，前m个元素留下来的概率均为：1/m * (1 - 1/(m+1)) = 1/(m+1)，也就是1/n。
6. 综上可知，算法成立。

 

问题二

将问题一中的条件变为，k为任意整数的情况，即要求最终返回的元素有k个，这就是水塘抽样（Reservoir Sampling）问题。要求是：取到第n个元素时，前n个元素被留下的几率相等，即k/n。

算法同上面思路类似，将1/n换乘k/n即可。在取第n个数据的时候，我们生成一个0到1的随机数p，如果p小于k/n，替换池中任意一个为第n个数。大于k/n，继续保留前面的数。直到数据流结束，返回此k个数。但是为了保证计算机计算分数额准确性，一般是生成一个0到n的随机数，跟k相比，道理是一样的。

同样采用归纳法来分析：

1. 初始情况 n <= k：此时每个元素留下的概率均为1。
2. 当 n = k+1 时，第k+1个元素留下的概率为k/(k+1)，前k个元素留下的概率均为：k/k * (1 - k/(k+1) * 1/k) = k/(k+1)，即上一步留下的概率乘以这一步留下的概率。
3. 假设 n = m 时，每个元素留下的概率均为 k/n = k/m。
4. 那么，当 n = m+1 时，第m+1个元素留下的概率为1/(m+1)，前m个元素留下的概率均为：k/m * (1 - k/(m+1) * 1/k) = k/(m+1)，其中：k/m为上一步留下来的概率，k/(m+1) * 1/k 为这一步不能留下来的概率（第m+1个留下来，同时池中一个元素被踢出的概率）。
5. 综上可知，算法成立。

 

伪代码如下：

//stream代表数据流
//reservoir代表返回长度为k的池塘

//从stream中取前k个放入reservoir；
for ( int i = 1; i < k; i++)
    reservoir[i] = stream[i];
for (i = k; stream != null; i++) {
    p = random(0, i);
    if (p < k) reservoir[p] = stream[i];
return reservoir;