How to efficiently choose $n$ subset out of a set of $m$ many numbers, in a randomized uniform manner?

Question

Problems:

It is fairly simple: we have a list of numbers $x_1, x_2, \ldots,x_n,\ldots, x_m$. Our goal is to randomly and uniformly choose a subset of $n$ many numbers out of these.

This means that, for any $i \in \{1,2,\ldots,n,\ldots,m\}$, the probability of choosing the value of $x_i$ must be: $$ n\frac{1}{m} = \frac{n}{m} $$

An algorithm:

An instruction is $\text{swap}(x_{\text{left}},x_{\text{right}})$. It simply swaps their values. I.e.

1. $t := x_{\text{left}}$
2. $x_{\text{left}} := x_{\text{right}}$
3. $x_{\text{right}} := t$

A suggested algorithm is: for each $i \in \{1,2,\ldots,n\}$, randomly and uniformly choose some $k_i \in \{1,2,\ldots,n, \ldots, m\}$, and then execute $\text{swap}(x_i, x_{k_i})$. Then return $x_1, x_2, \ldots, x_n$ as the set of $n$-many chosen numbers.

The challenge:

Intuitively, that algorithm looks to me to be perfectly fine. I see absolutely no problem in it.

But when I try to look at it mathematically, I fail to prove it. Instead, I actually seem to be prove that it is not a correct solution!

First, let's look at the probability of choosing the first $n$ numbers:

For any $i \in \{1,2,\ldots,n\}$, $x_i$ can be chosen if:

• $x_i$ exists in the right hand of $\text{swap}$. There are exactly $n$ cases, including the case when $x_i$ exists in, both, the left and the right hands.
• $x_i$ exists in the left hand, such that there exists $x_j$ in the right hand such that $j \in \{1,2,\ldots,n\}$. There are exactly $n$ such cases, one of which is the case when $i=j$ which we have counted earlier. Therefore, to avoid counting the case of $i=j$ twice, we assume that there are $n-1$ cases.

There is no other case where $x_i$ is chosen. Therefore the probability of choosing a number $x_i$, given that $i \in \{1,2,\ldots,n\}$ is: $$ \frac{n + (n-1)}{m^2} $$

Now, let's look at the probability of choosing the reset of the numbers up to $m$:

For any $i \in \{n+1,n+2,\ldots,m\}$, $x_i$ can be chosen if:

• $x_i$ exists in the right hand, such that there exists $x_j$ in the right hand such that $j \in \{1,2,\ldots,n\}$. There are exactly $n$ such cases.

There is no other case of having $x_i$ chosen. So the probability of choosing $x_i$ is: $$ \frac{n}{m^2} $$

Comparing the probabilities:

That says that there is a bias. The probability of choosing the 1st $n$ numbers is higher than the probability of choosing the reset of the numbers up to the $m^{th}$

My question:

Where did I go wrong? How to address this problem?

Answer 1

The "suggested algorithm" is incorrect. One way to see why that is so is to count the number of equiprobable permutations performed by the algorithm. At each step there are $m$ possible values for $k$, whence after $n$ steps there are $m^n$ possible results. Although many results will be duplicated, the point is that the probability of outputting any particular permutation must be a multiple of $m^{-n}$. However, the correct probability, $1/m!$, is rarely such a multiple. For instance, with $m=3$ and $n=2$ you will produce $3^2=9$ possible permutations, each with chance $1/9$, but there are only $m!=6$ distinct permutations. Since $1/6$ is not a multiple of $1/9$, none of the permutations will be produced with a chance of $1/6$.

There are better ways. One is "Algorithm P" from Knuth's Seminumerical Algorithms, section 3.4.2:

for i := 1 to n do
    Swap(x[i], x[RandInt(i,m)])

The proof that this works is by induction on $m$.

1. Obviously it works in the base case (either $m=0$ or $m=1$, as you prefer).

2. In the first step, every $x_k$ has an equal chance of being moved into the first position. The algorithm proceeds to work recursively and independently on the elements in positions $2$ through $n$, where by the inductive hypothesis every one of those elements has an equal and independent chance of appearing anywhere among the first $n-1$ positions. Consequently all $m$ of the elements of $x$ have equal and independent chances of appearing among the first $n$ positions when the algorithm terminates, QED.

Knuth's Algorithm S guarantees the output will be in the order in which they originally appeared in $x$:

Select := n
Remaining := m
for i := 1 to m do
    if RandReal(0,1) * Remaining < Select then
        output x[i]
        Select--
    Remaining--

Exercises (2), (3), and (4) in that section ask the reader to prove this algorithm works. Once again the proof is an induction on $m$.

1. When $m=n=1$, $x$ itself is always returned.

2. Otherwise, $x_1$ will be output with probability $n/m$ in the first step, which is the correct probability, and the algorithm proceeds recursively to output $n-1$ elements of $x_{-1} = (x_2, x_3, \ldots, x_m)$ in sorted order if $x_1$ was output and otherwise will output $n$ elements of $x_{-1}$ in sorted order, QED.

---

If you would like to follow along, here is an executable version in R, along with a quick simulation to verify that all elements of $x$ have equal chances of being included.

algorithm.S <- function(x, n=1) {
  m <- length(x)
  #
  # Check input.
  #
  if (n < 0 || n > m) stop("Subset size out of range.")
  #
  # Handle special cases that R has trouble with.
  #
  if (m <= 1) 
    if (n==0) return (x[c()]) else return(x)
  #
  # The algorithm.
  #
  y <- x[1:n]
  select <- n
  remaining <- m
  j <- 0
  for (i in 1:m) {
    if (runif(1) * remaining < select) {
      j <- j+1; y[j] <- i
      select <- select-1
    }
    remaining <- remaining-1
  }
  return(y)
}

x <- 1:10
sim <- replicate(1e4, algorithm.S(x, 4))
hist(sim, breaks=seq(min(x)-1/2, max(x)+1/2, by=1))

Answer 2

Start at the beginning of the list.

Pick the element with probability equal to $\frac{n}{m}$. If it is chosen set $n = n - 1$. Set $m = m - 1$ Now pass on to the next element and repeat until either $n$ or $m$ is zero.

Not tested in detail but should work.

Or use the facilities of your favourite statistical software.