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. 参考资料
蓄水池抽样——《编程珠玑》读书笔记