An alternative sampling without replacement

Question

Consider a set $X := \{x_1, \ldots, x_n\}$ with corresponding weights $p_1, \ldots, p_n$. Suppose we would like to draw $m < n$ distinct (i.e. unique) elements in a way that the probability of drawing a set $(x_{i_1}, \ldots, x_{i_m})$ is proportional to the product of corresponding weights, i.e. \begin{equation}\tag{1} p_{(i_1, \ldots , i_m)} :\propto p_{i_1} \cdots p_{i_m}. \end{equation} Since we are drawing among all sequences $(i_1, \ldots, i_m)$ with no duplicates, the normalizing constant is the sum over all such sequences: $$\tag{2} p_{(i_1, \ldots , i_m)} = \frac{p_{i_1} \cdots p_{i_m}}{\sum_{\substack{(j_1, \ldots, j_m) \\ \text{no duplicates}}} p_{j_1} \cdots p_{j_m}}. $$ This is equivalent to a sampling scheme where one samples $m$ weighted elements with replacement and then rejects (and repeats) the entire sample if it contains at least one pair duplicate values. Surprisingly (to me) this scheme is not equivalent to weighted sampling without replacement, i.e. where one sequentially samples first $i_1$ with weights, updates all remaining weights and then samples $i_2$, etc.

When $n$ is large, direct sampling by computing the denominator in (2) is practically impossible. Further, if $m$ is close to $n$ or if weights have a very non-uniform distribution, the rejection sampler may take very long before the sampling with replacement manages to draw a sample of $m$ distinct draws. Yet (1) has a simple structure, which perhaps could enable more efficient sampling.

My question is: Is there a way to sample efficiently from this scheme? This could for instance be a sequential sample with some update of weights.

Answer 1

There are algorithms in the survey literature for what is called "$\pi$ps sampling" or "PPS sampling without replacement". Most of these only get the marginal probabilities for individual observations to match what you want and the pairwise probabilities to be non-zero (which is still not entirely trivial). The Conditional Poisson Sampling approach gets the joint probabilities as you want them (I think). It uses rejection, but in a different way.

The idea is to take probabilities $q_i$ and independently sample each observation or not with probability $q_i$, then accept the sample if the sample size turns out to be $m$. The probability of getting the right sample size is going to be of order $m^{-1/2}$ -- if the probabilities were the same, the sample size would be Poisson with variance $m$, so the number of high-probability values it has would be of the order of $\sqrt{m}$.

The problem is computing the $q_i$. As you note, they are not the $p_i$, and computing them does seem to require computing the big messy denominator. It's not completely intractable if $m$ is small, but it is slow. This is one reference.

A good book about algorithms for PPS sampling is the one by Yves Tillé; he also has software