Algorithm for sampling fixed number of samples from a finite population

Question

I'm looking for an algorithm that would do the following:

Imagine that you need to sample uniformly at random and without replacement $k$ elements from a pool of $n$ elements. The catch is that $n$ is unknown. You can iterate over the whole set of the elements, so eventually, you would learn $n$. The algorithm is ought to have $O(n)$ time complexity, you are allowed to iterate over the set only once. You are allowed to keep a small cache of size $m \approx k$, but memorizing all the elements is not possible as $n$ could be very large.

Example: you work at a factory, your job is to select exactly $k$ items from the production line for an inspection per day. You need to sample them uniformly at random (with a probability $k/n$ each). The problem is that the number of items produced by the production line is unpredictable, so you don't know if it would produce $k$ items, or $1,000k$ items, or a completely different number. You need to decide on the fly to keep an item or leave it, as they move over the production line and won't be available later.

This should be possible. I guess what I can do is to take the first $k$ samples that I saw, and on each following step with some probability $p$ drop one of the already available samples and take instead a new one. The probability of drawing a new sample would need to evolve over time, but also the probabilities of rejecting the already sampled values would change. Probability of rejecting $x_t$ would need to depend on time order when it was observed $t$. The first $k$ values were sampled with probability $1$, $k+1$ value would be sampled with probability $k/k+1$, the last value with probability $k/n$, etc. By rejecting items already collected we would "correct" the sampling probabilities.

I don't want to re-invent the wheel, so I wanted to ask if there already is an algorithm like this?

Answer 1

Yes.

Collect the first $k$ items encountered into the cache. At steps $j=k+1, \ldots, n,$ place item $j$ in the cache with probability $k/j,$ in which case you will remove one of the existing items uniformly at random. After you have been through the entire population, the cache will be the desired random sample.

This algorithm is similar to a standard algorithm for creating a random permutation of $n$ items. It's essentially Durstenfeld's version of the Fisher-Yates shuffle.

Here is a diagram of how such a sample of size $k=20$ evolved for a population that eventually was size $n=300.$ The lines at each iteration indicate the indexes of the sample members.

At each iteration, the sample should be roughly uniformly distributed between $1$ and the iteration--conditional, of course, on how uniformly distributed it had previously been. Of crucial importance is to note how some of the earliest elements (shown in red) manage to persist in the sample to the end: these need to have the same chances of being in the sample as any of the later elements.

---

To prove the algorithm works, we may view it as a Markov chain.

• The set of states after $n\ge k$ items have been processed can be identified with the set of $k$-subsets $\mathcal{I} = \{i_1, i_2, \ldots, i_k\}$ of the indexes $1,2,\ldots, n$ denoting which items are currently in the sample.

• The algorithm makes a random transition from any subset $\mathcal I$ of $\{1,2,\ldots, n\}$ to $k+1$ distinct possible subsets of $\{1,2,\ldots, n, n+1\}.$ One of them is $\mathcal I$ itself, which occurs with probability $1 - k/(n+1).$ The other of them are the subsets where $i_j$ is replaced by $m+1$ for $j=1,2,\ldots, k.$ Each of these transitions occurs with probability

$$\frac{1}{k}\left(\frac{k}{n+1}\right) = \frac{1}{n+1}.$$

We need to prove that after $n \ge k$ steps, every $k$-subset of $\{1,2,\ldots, n\}$ has the same chance of being the sample. We can do this inductively. To this end, suppose after step $n\ge k$ that all $k$ subsets have equal chances of being the sample. These chances therefore are all $1/\binom{n}{k}.$ After step $n+1,$ a given subset $\mathcal I$ of $\{1,2,\ldots, n+1\}$ can have arisen as a transition from $n-k+2$ subsets of $\{1,2,\ldots, n\}:$ namely,

• If $\mathcal{I}$ does not contain $n+1,$ it arose as a transition of probability $1-k/(n+1)$ from itself, where it originally had a chance of $1/\binom{n}{k}$ of occurring. Such subsets therefore appear with individual chances of

  $$\Pr(\mathcal{I}) = \frac{1}{\binom{n}{k}} \times \left(1 - \frac{k}{n+1}\right) = \frac{1}{\binom{n+1}{k}}.$$

• If $\mathcal{I}$ does contain $n+1,$ it arose upon replacing one of the $n-(k-1)$ indexes in $\{1,2,\ldots, n\}$ that do not appear in $\mathcal I$ with the new index $n+1.$ Each such transition occurs with chance $1/(n+1),$ again giving a total chance of

  $$\Pr(\mathcal I) = (n-(k-1)) \times \frac{1}{\binom{n}{k}} \times \frac{1}{n+1} = \frac{1}{\binom{n+1}{k}}.$$

Consequently, all possible $k$-subsets of the first $n+1$ indexes have a common chance of $1/\binom{n+1}{k}$ of occurring, proving the induction step.

To start the induction, notice that at step $n=k$ there is exactly one subset and it has the correct chance of $1$ to be the sample! This completes the proof.

---

This R code demonstrates the practicality of the algorithm. In the actual application you would not have a full vector population: instead of looping over seq_along(population), you would have a source of data from which you sequentially fetch the next element (as in population[j]) and increment j until it is exhausted.

sample.online <- function(k, population) {
  cache <- rep(NA, k)
  for (j in seq_along(population)) {
    if (j <= k) {
      cache[j] <- population[j]
    } else {
      if (runif(1, 0, j) <= k) cache[sample.int(k, 1)] <- population[j]
    }
  }
  cache
}

Answer 2

I think the most intuitive solution is that you have an ordered list, and every time you see a new item, you place the item into the list at a random location. Then you take the first $k$ elements of that list.

Since you're taking the first $k$ elements, you don't need to keep track of the elements after that, so you can instead maintain a list of length $k$. Each time you see a new item, if you had seen $m$ items previously, there's $m+1$ different places to put the new item, so each for each slot in the partial list, there's a $\frac 1 {m+1}$ probability of the new item going there, and a $\frac{m+1-k}{m+1}$ probability of it not going in the partial list at all. If it does go in a partial list, then item that was previously last in the list gets pushed out.

Since the order doesn't matter, you can simplify it even further by having each new item having probability $\frac k {m+1}$ of being added, and if it is added, then each old item has probability $ \frac 1 k$ being dropped.

This is the algorithm in whuber's answer (my $m+1$ is their $j$), but I think the explanation for it is more intuitive.

Answer 3

Complementary to @whuber's and @Accumulation's good answers (+1).

Sampling techniques used to address such tasks are usually categorised under the umbrella of reservoir sampling; these sampling methodologies have been strongly motivated by the need to sample streaming data where the overall sample size $n$ is unknown or by definition dynamic. The term "reservoir" refers to the size of the resulting sample. The original principle is discussed in "Random sampling with a reservoir" (1985) by Vitter but came to prominence with social media applications; "TeRec: a temporal recommender system over tweet stream" (2013) by Chen et al. is a short and straightforward implementation of how the original algorithm of Vitter was adapted/extended to suit the needs of social media apps like Twitter and Weibo. I came also came across an excellent blog post on the matter here: Reservoir sampling by Startin which takes a much more programmatic approach.

Please note that the original request from Tim asked for $O(n)$ complexity but they are algorithms that can do even better than that. Similarly, they are implementations that offer the ability to do weighted reservoir sampling. (i.e. the probability of each item to be selected is determined by its relative weight - see "Weighted random sampling with a reservoir" (2005) by Efraimidis and Spirakis for an early work on that)