21
$\begingroup$

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?

$\endgroup$
1
  • 7
    $\begingroup$ I might be missing something here but isn't this a simple application of reservoir sampling? And you can go below $O(n)$ if needed. $\endgroup$
    – usεr11852
    Mar 30, 2022 at 1:02

3 Answers 3

23
$\begingroup$

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.

Figure

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
}
$\endgroup$
8
  • 2
    $\begingroup$ Seems like I was overthinking it. It's beautiful how simple it is. $\endgroup$
    – Tim
    Mar 29, 2022 at 22:20
  • 8
    $\begingroup$ What caught my interest is that it is adaptive: one could begin with a SRS from a population, later discover that the entire population was not available at the time it was originally sampled, and use this algorithm to modify the original sample into an SRS (of the same size) from the increased population. This might have saved me considerable trouble last year when sampling medical records in a legal case: the parties kept coming up with more records as the case progressed. The analyses I did were accordingly very complicated and an approach like this one could have simplified it. $\endgroup$
    – whuber
    Mar 29, 2022 at 22:24
  • 3
    $\begingroup$ (+1) As I say my comment on the main question, there is a substantial body of work on this matter under reservoir sampling (as I interpret this application at least). It has been strongly motivated by the need to sample streaming data. I would discourage anyone from implementing such an algorithm from scratch aside from pedagogical reasons; if for example, we are interested in R, there is the package stream that seems to do this work for us. Other languages will have similar implementation. $\endgroup$
    – usεr11852
    Mar 30, 2022 at 8:56
  • 3
    $\begingroup$ @usεr11852 Thank you for those observations! I was unaware of any of that work. I proposed this algorithm directly in response to the question, relying only on basic knowledge of algorithms and data structures. I am pleased now to know this is "Algorithm R" and that its natural generalizations (skipping and weighting) have been published. (I do want to observe that the Wikipedia article is incomplete because it only demonstrates the sample is uniform.) I don't see any problems with coding this algorithm from scratch, though: it is simple and easily checked with simulations (which I did). $\endgroup$
    – whuber
    Mar 30, 2022 at 14:24
  • 3
    $\begingroup$ @Tim: I partially agree. I think whuber and yourself are not "average blokes" when it comes to your understanding of implementing a statistical procedure. If anything, he is an expert and even wrote a unit test for his own code. For that matter, if I had to put a piece of code in production/report I would prefer something that has had more eyeballs on it than just my own (or any other "single pair of eyeballs"). $\endgroup$
    – usεr11852
    Mar 30, 2022 at 21:31
10
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ +1 nice intuition, plus it reveals possible code optimization: direct implementation of whuber’s algorithm needs a coin toss (pick value or ignore) and $k$-sided dice roll (where to insert it), while it can be done with $m+1$ sided dice alone. $\endgroup$
    – Tim
    Mar 30, 2022 at 7:46
  • $\begingroup$ +1 I like this explanation. $\endgroup$
    – whuber
    Mar 30, 2022 at 22:25
9
$\begingroup$

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)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.