Sampling a sequence without replacement with non-uniform probability of sampling each element Suppose we have N objects which we wish to sample from, and each object $i$ has a certain probability of being chosen, where of course $\sum_{i=1}^N p_i = 1$.
We wish to sample without replacement.
Conditioned on object $i$ being chosen first (with probability $p_i$), the probability of choosing object $j$ second will be the renormalized probabilty
$$
\frac{p_j}{\sum_{l\neq i} p_l}\,.
$$
The unconditional probability of object $j$ being chosen second will then be
$$
\sum_{i\neq j} p_i \frac{p_j}{\sum_{l\neq i} p_l}\,.
$$
The unconditional probability of object $j$ being chosen third will then be
$$
\sum_{i,m \ s.t. \ i\neq m\neq j} p_i \frac{p_m}{\sum_{l\neq i} p_l} \frac{p_j}{\sum_{n\notin \{i, m\}} p_n} \,.
$$
The probability of object $j$ being chosen on the $k$th trial is then in general a pretty complicated expression...
Is there a way to simplify this to something "nice", or can we make some approximation (maybe for large N)? Or can we at the very least come up with an expression for the expected number of trials before we choose object $j$? It feels like we could ~almost~ use the negative hypergeometric distribution, but not quite... Any ideas?
 A: There is a very clever way to sample a probability distribution without replacement, without having to renormalize.
Recall that the exponential distribution with parameter $\lambda$ is the random variable with CDF
$$
F_{\lambda}(x)=1-e^{-\lambda x},\qquad x\ge 0
$$
It is easy to generate a random exponential sample given uniformly random bits, using inverse transform sampling.
Now, to draw a sample of size $s$ from a population of size $N$ with weights $p_1,\dots,p_N$, do the following.

*

*Generate $N$ independent random numbers, $X_1,X_2,\dots,X_N$, where $X_i$ is exponential with parameter $p_i$.


*Sort these $N$ numbers from smallest to largest, remembering their indices.


*The indices of the $s$ smallest numbers correspond to the selected items.
This procedure takes $O(N\log N)$ time. By contrast, naive renormalization requires at least $\Omega(N^2)$ time. The tradeoff is that this procedure is technically an approximation. If the random variables $X_i$ were stored with infinite precision, this would be exactly correct, but quantization errors will make the probabilities a little different.

To prove that this procedure accurately simulates the "weighted-sampling-without-replacement" procedure, you need only these three facts about the exponential distribution.

*

*If $X$ is exponential with parameter $\lambda$, and $Y$ is independently exponential with parameter $\mu$, then $\min(X,Y)$ is exponential with parameter $\lambda +\mu$.


*With the same assumptions as the last point, $P(X<Y)=\frac{\lambda}{\lambda +\mu}$.


*The exponential distribution is memoryless, in the sense that $P(X\le s+t\mid X>t)=P(X\le s)$ for all $s,t\ge 0$.
This implies that, for example, the probability that the first sampled element is $X_1$ is $P(X_1<\min(X_2,X_3,\dots,X_N))$, which is $\frac{p_1}{p_1+(p_2+\dots+p_N)}$, and then for subsequent samples, you can do a similar calculation where you leave $p_1$ out (because of memorylessness).
