Reservoir sampling algorithm probability

Question

I'm reading about the reservoir sampling technique called Algorithm R.

The idea is we can take a sample of size $n$ from a population of size $N$ even when $N$ is unknown/too expensive to retrieve in $N$ time. I quote a sample implementation from wikipedia:

(*
 S has items to sample, R will contain the result
*)
ReservoirSample(S[1..n], R[1..k])
 // fill the reservoir array
 for i = 1 to k
     R[i] := S[i]

 // replace elements with gradually decreasing probability
 for i = k+1 to n
   j := random(1, i)   // important: inclusive range
   if j <= k
       R[j] := S[i]

The explanation for why this works:

The algorithm creates a "reservoir" array of size $k$ and populates it with the first $k$ items of $S$.

That is pretty clear, the first for loop in the code sample does that. Then:

It then iterates through the remaining elements of $S$ until $S$ is exhausted. At the $i$-th element of $S$, the algorithm generates a random number $j$ between $1$ and $i$. If $j$ is less than or equal to $k$, the $j$-th element of the reservoir array is replaced with the $i$-th element of $S$.

Also clear. But here comes the Math:

In effect, for all $i$, the $i$-th element of $S$ is chosen to be included in the reservoir with probability $\frac{k}{i}$.

Hmm, I'm probably not well versed enougn in probablility theory, so can someone explain in more detail how the probability is $\frac{k}{i}$?

What I also don't get: do all elements have the same chance of being included in the (final) reservoir? From my understanding, it seems like the probability decreases for the elements at the end of $S$ to be included, then the reservoir will be skewed towards including the earlier elements, or am I misled?

Answer 1

The whole reason for performing this sampling method is to get an uniform sample even if the population size is unknown at the start. So, if this method works, the probability cannot be skewed.

What happens, is that the samples chosen in the beginning have a chance to be overwritten by the samples chosen later. Intuitively, this balances it somewhat. But we must show that this 'balance' holds precisely:

The probability of $\frac{k}{i}$ refers to the probability that sample $i$ is added to the reservoir, as it is added when a uniform random integer from $[1,i]$ is less than $k$, i.e. $k$ positions in the range $[1,i]$ lead to inclusion in the reservoir. (Total number of events is $i$, number of events that we 'want' is $k$, all events are equally likely, so final probability is $\frac{k}{i}$)

Now what is the probability that $i$, when chosen, remains in the reservoir until the end? It gets removed when a single integer in the range $[1,m]$ is picked, for all $i<m\leq N$. The probability that it doesn't get replaced at position $m$ is $1-\frac{1}{m}$. As these events are independent, the total probability that $i$ doesn't get replaced is $\prod_{m=i+1}^N (1-\frac{1}{m}) = \frac{i}{N}$ $^*$. As $i$ is in the reservoir at the end if it gets picked and not replaced, this probability is $\frac{k}{i}\cdot \frac{i}{N} = \frac{k}{N}$.

We're not done yet, as the calculation above only holds for the samples not originally in the reservoir. However, that case is similar: as they are always put in the reservoir, we only need the probability they stay, which is $\prod_{m=k+1}^N (1-\frac{1}{m}) = \frac{k}{N}$.

So the probability of being included in the final sample is the same for all positions, which means the final sample is uniform.

*: The expression $\prod_{m=i+1}^N (1-\frac{1}{m})$ can be simplified as follows:

\begin{align}\prod_{m=i+1}^N (1-\frac{1}{m}) &= \prod_{m=i+1}^N (\frac{m-1}{m})\\ &= \frac{i}{i+1}\cdot\frac{i+1}{i+2}\cdots \frac{N-2}{N-1}\cdot \frac{N-1}{N}\\ &= i\cdot\frac{i+1}{i+1}\cdot\frac{i+2}{i+2} \cdots \frac{N-1}{N-1}\cdot \frac{1}{N} =\frac{i}{N}. \end{align}

It is also possible to show by induction on $N$, but I think the previous is more insightful.