蓄水池取样(转)

        问题:如何随机从n个对象中选择一个对象,这n个对象是按序排列的,但是在此之前你是不知道n的值的。

        思路:如果我们知道n的值,那么问题就可以简单的用一个大随机数rand()%n得到一个确切的随机位置,那么该位置的对象就是所求的对象,选中的概率是1/n。

        但现在我们并不知道n的值,这个问题便抽象为蓄水池抽样问题,即从一个包含n个对象的列表S中随机选取k个对象,n为一个非常大或者不知道的值。通常情况下,n是一个非常大的值,大到无法一次性把所有列表S中的对象都放到内存中。我们这个问题是蓄水池抽样问题的一个特例,即k=1。

        解法:我们总是选择第一个对象,以1/2的概率选择第二个,以1/3的概率选择第三个,以此类推,以1/m的概率选择第m个对象。当该过程结束时,每一个对象具有相同的选中概率,即1/n,证明如下。

        证明:第m个对象最终被选中的概率P=选择m的概率*其后面所有对象不被选择的概率,即

        对应的该问题的伪代码如下:

  1. i = 0
  2. while more input items
  3. with probability 1.0 / ++i
  4. choice = this input item
  5. print choice

        C++代码实现如下:

  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <ctime>
  4. #include <vector>
  5. using namespace std;
  6. typedef vector<int> IntVec;
  7. typedef typename IntVec::iterator Iter;
  8. typedef typename IntVec::const_iterator Const_Iter;
  9. // generate a random number between i and k,
  10. // both i and k are inclusive.
  11. int randint(int i, int k)
  12. {
  13. if (i > k)
  14. {
  15. int t = i; i = k; k = t; // swap
  16. }
  17. int ret = i + rand() % (k - i + 1);
  18. return ret;
  19. }
  20. // take 1 sample to result from input of unknown n items.
  21. bool reservoir_sampling(const IntVec &input, int &result)
  22. {
  23. srand(time(NULL));
  24. if (input.size() <= 0)
  25. return false;
  26. Const_Iter iter = input.begin();
  27. result = *iter++;
  28. for (int i = 1; iter != input.end(); ++iter, ++i)
  29. {
  30. int j = randint(0, i);
  31. if (j < 1)
  32. result = *iter;
  33. }
  34. return true;
  35. }
  36. int main()
  37. {
  38. const int n = 10;
  39. IntVec input(n);
  40. int result = 0;
  41. for (int i = 0; i != n; ++i)
  42. input[i] = i;
  43. if (reservoir_sampling(input, result))
  44. cout << result << endl;
  45. return 0;
  46. }

        对应蓄水池抽样问题,可以类似的思路解决。先把读到的前k个对象放入“水库”,对于第k+1个对象开始,以k/(k+1)的概率选择该对象,以k/(k+2)的概率选择第k+2个对象,以此类推,以k/m的概率选择第m个对象(m>k)。如果m被选中,则随机替换水库中的一个对象。最终每个对象被选中的概率均为k/n,证明如下。

        证明:第m个对象被选中的概率=选择m的概率*(其后元素不被选择的概率+其后元素被选择的概率*不替换第m个对象的概率),即


        蓄水池抽样问题的伪代码如下:

  1. array S[n]; //source, 0-based
  2. array R[k]; // result, 0-based
  3. integer i, j;
  4. // fill the reservoir array
  5. for each i in 0 to k - 1 do
  6. R[i] = S[i]
  7. done;
  8. // replace elements with gradually decreasing probability
  9. for each i in k to n do
  10. j = random(0, i); // important: inclusive range
  11. if j < k then
  12. R[j] = S[i]
  13. fi
  14. done

        C++代码实现如下,该版本假设n知道,但n非常大:

  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <ctime>
  4. using namespace std;
  5. // generate a random number between i and k,
  6. // both i and k are inclusive.
  7. int randint(int i, int k)
  8. {
  9. if (i > k)
  10. {
  11. int t = i; i = k; k = t; // swap
  12. }
  13. int ret = i + rand() % (k - i + 1);
  14. return ret;
  15. }
  16. // take m samples to result from input of n items.
  17. bool reservoir_sampling(const int *input, int n, int *result, int m)
  18. {
  19. srand(time(NULL));
  20. if (n < m || input == NULL || result == NULL)
  21. return false;
  22. for (int i = 0; i != m; ++i)
  23. result[i] = input[i];
  24. for (int i = m; i != n; ++i)
  25. {
  26. int j = randint(0, i);
  27. if (j < m)
  28. result[j] = input[i];
  29. }
  30. return true;
  31. }
  32. int main()
  33. {
  34. const int n = 100;
  35. const int m = 10;
  36. int input[n];
  37. int result[m];
  38. for (int i = 0; i != n; ++i)
  39. input[i] = i;
  40. if (reservoir_sampling(input, n, result, m))
  41. for (int i = 0; i != m; ++i)
  42. cout << result[i] << " ";
  43. cout << endl;
  44. return 0;
  45. }

        下面这个程序假设不知道n的大小

  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <ctime>
  4. #include <vector>
  5. using namespace std;
  6. typedef vector<int> IntVec;
  7. typedef typename IntVec::iterator Iter;
  8. typedef typename IntVec::const_iterator Const_Iter;
  9. // generate a random number between i and k,
  10. // both i and k are inclusive.
  11. int randint(int i, int k)
  12. {
  13. if (i > k)
  14. {
  15. int t = i; i = k; k = t; // swap
  16. }
  17. int ret = i + rand() % (k - i + 1);
  18. return ret;
  19. }
  20. // take m samples to result from input of n items.
  21. bool reservoir_sampling(const IntVec &input, IntVec &result, int m)
  22. {
  23. srand(time(NULL));
  24. if (input.size() < m)
  25. return false;
  26. result.resize(m);
  27. Const_Iter iter = input.begin();
  28. for (int i = 0; i != m; ++i)
  29. result[i] = *iter++;
  30. for (int i = m; iter != input.end(); ++i, ++iter)
  31. {
  32. int j = randint(0, i);
  33. if (j < m)
  34. result[j] = *iter;
  35. }
  36. return true;
  37. }
  38. int main()
  39. {
  40. const int n = 100;
  41. const int m = 10;
  42. IntVec input(n), result(m);
  43. for (int i = 0; i != n; ++i)
  44. input[i] = i;
  45. if (reservoir_sampling(input, result, m))
  46. for (int i = 0; i != m; ++i)
  47. cout << result[i] << " ";
  48. cout << endl;
  49. return 0;
  50. }

        本文参考:

http://www.cnblogs.com/HappyAngel/archive/2011/02/07/1949762.html

http://en.wikipedia.org/wiki/Reservoir_sampling

http://en.wikipedia.org/wiki/Fisher-Yates_shuffle


原文地址:https://www.cnblogs.com/developing/p/11160928.html