C++: Knuth’s Algorithm S

This is a method of randomly sampling n items from a set of M items, with equal probability; where M >= n and M, the number of items is unknown until the end. This means that the equal probability sampling should be maintained for all successive items > n as they become available (although the content of successive samples can change).

The algorithm

1. Select the first n items as the sample as they become available;
2. For the i-th item where i > n, have a random chance of n/i of keeping it. If failing this chance, the sample remains the same. If not, have it randomly (1/n) replace one of the previously selected n items of the sample.
3. Repeat #2 for any subsequent items.

The Task

1. Create a function s_of_n_creator that given n the maximum sample size, returns a function s_of_n that takes one parameter, item.
2. Function s_of_n when called with successive items returns an equi-weighted random sample of up to n of its items so far, each time it is called, calculated using Knuths Algorithm S.
3. Test your functions by printing and showing the frequency of occurrences of the selected digits from 100,000 repetitions of:

1. Use the s_of_n_creator with n == 3 to generate an s_of_n.
2. call s_of_n with each of the digits 0 to 9 in order, keeping the returned three digits of its random sampling from its last call with argument item=9.

Note: A class taking n and generating a callable instance/function might also be used.

Reference

• The Art of Computer Programming, Vol 2, 3.4.2 p.142

#include <iostream>
#include <functional>
#include <vector>
#include <cstdlib>
#include <ctime>
 
template <typename T>
std::function<std::vector<T>(T)> s_of_n_creator(int n) {
  std::vector<T> sample;
  int i = 0;
  return [=](T item) mutable {
    i++;
    if (i <= n) {
      sample.push_back(item);
    } else if (std::rand() % i < n) {
      sample[std::rand() % n] = item;
    }
    return sample;
  };
}
 
int main() {
  std::srand(std::time(NULL));
  int bin[10] = {0};
  for (int trial = 0; trial < 100000; trial++) {
    auto s_of_n = s_of_n_creator<int>(3);
    std::vector<int> sample;
    for (int i = 0; i < 10; i++)
      sample = s_of_n(i);
    for (int s : sample)
      bin[s]++;
  }
  for (int x : bin)
    std::cout << x << std::endl;
  return 0;
}

30052
29740
30197
30223
29857
29688
30095
29803
30098
30247

SOURCE

Content is available under GNU Free Documentation License 1.2.