Algorithm for sampling fixed number of samples from a finite population 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?

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
}

