# Sampling from probability distribution via multiplying probabilities by a uniform random number

When reading about a certain piece of software (NUS ScheduleTool, for generating nonuniform sampling schedules for biomolecular NMR spectroscopy), I found the following curious method for drawing $N$ unique points given a discrete probability distribution:

Suppose we have a discrete set $S$ with a probability mass function $p\colon S \to [0,1]$ (summing to $1$). For each point $x \in S$, multiply its probability mass $p(x)$ by a random number drawn from a uniform distribution on the unit interval $[0,1]$ to get a new number $g(x)$. Sort the points in descending order by their values of $g$ and pick the $N$ first ones (where $N \leq |S|$).

Does this process have anything to do with the usual method of sampling without replacement? Is it possible to relate the distribution of the drawn points to the original distribution in any simple way?

Comparing this to the numpy.random.choice implementation of conventional sampling without replacement, I found that the resulting point distributions do look noticeably different (in the case of an exponential distribution where $S$ is a 2D grid, anyway). So it seems that the methods really give different results (unless there is some issue with the way I'm using their software).

I do not know the software authors' motivation for using this method. Since they did not provide any reasons for it, I'm inclined to suspect that this was an ad-hoc attempt to reproduce conventional sampling, but a more charitable interpretation is that there is some good reason for doing it this way.

To clarify what I mean by conventional sampling without replacement: I assume that we know how to draw one item from a nonuniform distribution $p$ (for example, sending a uniform variable through the inverse CDF to get an element index). Then, sampling multiple items without replacement means:

• Draw the first point $x_1$ from $p_1 := p\colon S \to [0,1]$.
• Draw the second point $x_2$ from $p_2 \colon S \setminus \{x_1\} \to [0,1]$.
• Draw the third point $x_3$ from $p_3 \colon S \setminus \{x_1, x_2\} \to [0,1]$.
• ...

where, at each stage, we get a new distribution by removing the sample we just drew, and adjusting the normalization: $p_{i+1}(x) = \frac{p_i(x)}{1 - p_i(x_i)}$.

I'm pretty sure that this is equivalent to repeatedly sampling with replacement, discarding any duplicates: E.g. for the first step, if we've drawn $x_1$ we keep drawing until we get an $x_2 \neq x_1$. So for $x \neq x_1$: $$p_2(x) = p(x) + p(x_1) p(x) + p(x_1)^2 p(x) + \cdots = \frac{p(x)}{1 - p(x_1)},$$ the same as before.

To be fair, I'm not sure how "conventional" these methods are.

• You have values: 1 and 2 with equal probabilities. The result of multiplication is 0.5 and 1. What do you do next? How to select the 50 smallest values? – zoli Aug 10 '18 at 16:36
• @zoli To select 50 items you would need to start with at least 50 items. You only select each item at most once. – Elias Riedel Gårding Aug 10 '18 at 16:41
• What is a point $x$ then? – zoli Aug 10 '18 at 16:44
• @zoli I've made some clarifications, let me know if it's still unclear! – Elias Riedel Gårding Aug 10 '18 at 16:53
• @joriki You're right, especially since they constrain those "unnormalized" weights to lie between 0 and 1. – Elias Riedel Gårding Aug 10 '18 at 18:48

## 1 Answer

Imagine your set $S$ had just two points, $x$ and $y$, where $p(x)=\frac13$ and $p(y)=\frac23$. Now, suppose you used this method to draw a sample of size $1$. What is the probability the sample is $\{x\}$?

You would expect it to be $1/3$, right? However, this method would pick $\{x\}$ with probability $1/4$. Indeed, if one variable $U$ is uniform on $[0,\frac13]$, and another $V$ is uniform on $[0,\frac23]$, then $(U,V)$ is uniform on the rectangle $[0,\frac13]\times[0,\frac23]$, and the region where $U\ge V$ is a triangle of area $\frac12(\frac13)^2=\frac1{18}$, which is one fourth of the total area $\frac13\cdot\frac23=\frac29$.

I would say their sampling method is a heuristic, and a pretty flawed one at that.

Furthermore, I think there is only one random process consistent with the description "sample from a distribution without replacement." Generate a list $X_1,X_2\dots$ of independent samples from that distribution. Scan the list from beginning onwards, putting every value you see in a bucket unless an earlier equal value was already put in the bucket, until the bucket has $N$ items.

A good question would be if there is a clever, "one-shot" way to simulate a sample according to this process...

Edit: Answering my own question, there is such a method! For each $x\in S$, let $Z_x$ be an exponential random variable with rate $p(x)$. If you choose a subset of $S$ corresponding to the indices of the smallest $N$ values of $\{Z_x\mid x\in S\}$, then the resulting set exactly has the distribution described.

• Wow, you just beat me to it. That's pretty incredible for a two-hour-old question. – Brian Tung Aug 10 '18 at 17:18