I want to sample from the integers $\{1, \dots, k\}$ with probabilities $\{ p_i \}_{i=1}^k$, with replacement, until I see $m$ distinct elements (call that $n$ times).

You can view the distribution I want to sample from as a multinomial distribution, but instead of fixing the number of samples $n$, let $n$ be the first number that gives an exactly $m$-sparse vector. (Alternatively, the maximal number that gives an $m$-sparse vector would also be okay.)

Using an alias method, implementing sampling literally as written above takes $O(n + k)$ time.

So, whether this is efficient to run depends on how big $n$ is. If the probabilities are relatively uniform and/or $k \gg m$, this shouldn't be too bad. But in a pathological case where $k = m + 1$ and two of the $p_i$s are extremely small, this will take a very long time.

This is a version of the coupon collector's problem, of which I've found non-uniform variants, but I haven't found one where you only need to get $m < k$ of the varieties of coupons – and anyway, I don't know if that line of reasoning will help with a more efficient sampling algorithm.


  • Can you find, or bound, the distribution of $n$? Particularly in terms of $\max p_i / \min p_i$ or similar?
  • Is there a more efficient algorithm to perform this sampling?
  • 1
    $\begingroup$ This earlier question was about I think the same setting, though the currently accepted answer there assumes the $p_i$ are uniform. $\endgroup$
    – Danica
    May 7, 2017 at 12:11

1 Answer 1


TLDR just generalize the coupon collector techniques.

Suppose you have a discrete state space Markov chain that evolves on the space of all subsets of $\{1,\ldots,k\}$ that have size between $0$ and $m$ (inclusive). This has size $\binom{k}{0} + \binom{k}{1} + \binom{k}{2} + \cdots \binom{k}{m}$. In particular notice that when $k=m$, this is $2^k$.

Say time starts at $0$. $X_0$ is the empty set with probability $1$. The marginal/transition distribution of $P(X_1 = \cdot \mid X_0 = 0)$ is uniform over the $k$ singletons. $P(X_2 = \{j,k\} \mid X_1 = \{j\}) = p_k$ where $k\neq j$ and $P(X_2 = \{j\} \mid X_1 = \{j\}) = p_j$. If you write out the big ugly transition matrix, you'll see each row only has $k$ nonzero elements, because you can only do $k$ things at each time step.

Using that big ugly transition matrix, you might figure the transition matrix for $|X_t|$ (the cardinality/size of $X_t$). With this you can describe the stopping time of interest:

$$ n = \inf\{t : |X_t| = m\}. $$ Notice that $n \in \{m, m+1, \ldots\}$ and $$ P(n = j) = P(|X_n| = j \mid |X_{n-1}| = j-1) P(|X_{n-1}| = j-1) . $$ Both these factors could be coded up.

Regarding sampling, it's faster (but more memory-intensive) to sample $X_t$ or $|X_t|$ instead of the whole multinomial enchilada. The hard part is instantiating and storing the transition matrix, but it’s very straightforward (more so for $X_t$).

  • $\begingroup$ In terms of the last line: I see that coding this up would help determine the distribution of $n$, but I don't see how it would help with actual sampling versus just running the plain sampling. (In any case, the practical problem this question was for is no longer relevant to me seven years later, but interesting to think about again nonetheless.) $\endgroup$
    – Danica
    Apr 20, 2023 at 18:37
  • $\begingroup$ Yeah I saw the date and thought twice about answering this. RE: sampling, it depends on $k$ probably. Both algos sample a discrete rv, but my suggestion avoids another loop to check whether the new sample is unique. Small potatoes (if they even are potatoes) $\endgroup$
    – Taylor
    Apr 20, 2023 at 18:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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