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

## 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?

• Your question is how to do it or is it about validity of your algorithm? The simple solution is to sample m numbers from continuous uniform, then pick the n greatest values and use their positions to take the appropriate cases from your data. – Tim Aug 4 '16 at 11:32
• Both, but primarily how to do it in the simplest way. Your solution is interesting, I like thinking about it. So first I sample $m$ unique numbers. Then I search the $m$ numbers to identify the top-$n$ largest numbers (I think this implies sorting $m$ numbers). But before I continue to using their positions, isn't it already computationally more complicated as it requires me to essentially also sort $m$ many numbers? – caveman Aug 5 '16 at 2:43
• It's certainly not the most efficient algorithm but it can be implemented as one-liner in most high-level languages. – Tim Aug 5 '16 at 5:48

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))


• Thank you. Can you check this please bpaste.net/show/d799b1973bd6 --- it's the output for my "incorrect" algorithm, it shows $4$ numbers chosen from eight numbers $0,1,2,3,4,5,6,7$, along with their probabilities (so it's a histogram). I got this after generating permutations for $100,000,000$ times. Looks uniform to my eyes. Do you agree that this is also uniform? – caveman Aug 5 '16 at 3:08
• Never mind with my comment above. I finally understood why my algorithm is wrong (i.e. why using $i$ as the lower bound for the prng is needed). Thank you for your explanation, and your amazing proofs. Highly appreciate it! – caveman Aug 5 '16 at 11:14

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.

• So this will require an additional space where I store, or take note of, the chosen numbers. Any idea how will this be effectively different than my algorithm? What I do is that, instead of me choosing each element, I simply swap the chosen randomly. This way no extra space is needed to mark the chosen (but of course it is destructive as I change the position of numbers in the original list/array). – caveman Aug 5 '16 at 2:46