5
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Tim Aug 4 '16 at 11:32
  • $\begingroup$ 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? $\endgroup$ – caveman Aug 5 '16 at 2:43
  • 1
    $\begingroup$ It's certainly not the most efficient algorithm but it can be implemented as one-liner in most high-level languages. $\endgroup$ – Tim Aug 5 '16 at 5:48
6
$\begingroup$

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

enter image description here

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – caveman Aug 5 '16 at 3:08
  • 1
    $\begingroup$ 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! $\endgroup$ – caveman Aug 5 '16 at 11:14
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – caveman Aug 5 '16 at 2:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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