6. 蓄水池抽样

问题描述 ：随机从一个数据流中选取1个或k个数，保证每个数被选中的概率是相同的。数据流的长度 \(n\) 未知或者是非常大。

6.1. 随机选择1个数

在数据流中，依次以概率 \(1\) 选择第一个数，以概率 \(\frac{1}{2}\) 选择第二个数（替换已选中的数），…，依此类推，以概率 \(\frac{1}{m}\) 选择第 m 个数（替换已选中的数）。 结束时（遍历完了整个数据流），每个数被选中的概率都是 \(\frac{1}{n}\) 。证明:

第 m 个对象最终被选中的概率 = 选择第 m 个数的概率 x 后续所有数都不被选择的概率

即

\[P = \frac{1}{m} \times \left( \frac{m}{m+1} \times \frac{m+1}{m+2} \times \cdots \times \frac{n-1}{n} \right) = \frac{1}{n}.\]
 1#include <iostream>
 2#include <vector>
 3#include <utility> // swap
 4#include <ctime>
 5#include <cstdlib> // rand, srand
 6using namespace std;
 7
 8typedef vector<int> VecInt;
 9typedef VecInt::iterator Itr;
10typedef VecInt::const_iterator CItr;
11
12// 等概率产生区间 [a, b] 之间的随机数
13int randInt(int a, int b)
14{
15  if (a > b) swap(a, b);
16  return a + rand() % (b - a + 1);
17}
18
19bool sample(const VecInt data, int &result)
20{
21  if (data.size() <= 0) return false;
22
23  //srand(time(nullptr)); // 设置随机seed
24
25  CItr it = data.begin();
26  result = *it;
27  int m;
28  for (m = 1, it = data.begin() + 1; it != data.end(); ++m, ++it)
29  {
30    int ri = randInt(0, m); // ri < 1 的概率为 1/(m+1)
31    if (ri < 1) result = *it;
32  }
33  return true;
34}

6.2. 随机选择k个数

在数据流中，先把读到的前 k 个数放入“池”中，然后依次以概率 \(\frac{k}{k+1}\) 选择第 k+1 个数，以概率 \(\frac{k}{k+2}\) 选择第 k+2 个数，…， 以概率 \(\frac{k}{m}\) 选择第 m 个数（m > k）。如果某个数被选中，则 随机替换 “池”中的一个数。最终每个数被选中的概率都为 \(\frac{k}{n}\) 。 证明:

第 m 个对象最终被选中的概率 = 选择第 m 个数的概率 x（其后元素不被选择的概率 + 其后元素被选择的概率 x 不替换第 m 个数的概率）

即

\[\begin{split}P &=\ \frac{k}{m} \times \left[ \left( (1-\frac{k}{m+1}) + \frac{k}{m+1} \times \frac{k-1}{k} \right) \times \left( (1-\frac{k}{m+2}) + \frac{k}{m+2} \times \frac{k-1}{k} \right) \times \right. \\ & \ \quad \left. \cdots \times \left( (1-\frac{k}{n}) + \frac{k}{n} \times \frac{k-1}{k} \right) \right] \\ &=\ \frac{k}{m} \times \frac{m}{m+1} \times \frac{m+1}{m+2} \times \cdots \times \frac{n-1}{n} \\ &=\ \frac{k}{n}.\end{split}\]
 1#include <iostream>
 2#include <vector>
 3#include <utility> // swap
 4#include <ctime>
 5#include <cstdlib> // rand, srand
 6using namespace std;
 7
 8typedef vector<int> VecInt;
 9typedef VecInt::iterator Itr;
10typedef VecInt::const_iterator CItr;
11
12const int k = 10;
13int result[k];
14
15// 等概率产生区间 [a, b] 之间的随机数
16int randInt(int a, int b)
17{
18  if (a > b) swap(a, b);
19  return a + rand() % (b - a + 1);
20}
21
22bool sample(const VecInt data)
23{
24  if (data.size() < k) return false;
25
26  //srand(time(nullptr)); // 设置随机seed
27
28  CItr it = data.begin();
29  for(int m = 0; m < k; ++m) result[m] = *it++;
30
31  for (int m = k; it != data.end(); ++m, ++it)
32  {
33    int ri = randInt(0, m);
34    if (ri < k) result[ri] = *it; // ri < k 的概率为 k/(m+1)
35  }
36  return true;
37}

6.3. 参考资料

  1. 蓄水池抽样——《编程珠玑》读书笔记