Sampling from categorical distribution when the sorting step is ignored

Question

The commonly used algorithm to sample from categorical distribution is

1. sort the vector of cumulative probabilities $\boldsymbol{p} = (p_1,p_2,\dots,p_n)$ in decreasing order from the categories appearing with greatest, to the ones appearing with lowest probabilities, to obtain the permuted vector $\boldsymbol{p}_{\pi}$,
2. sample $u$ from uniform distribution on $(0,1)$,
3. accept the lowest $k$ such that $p_{\pi(k)} \ge u$ as $X$ (where $\pi(k)$ is position of $k$ in permuted vector $\boldsymbol{p}_{\pi}$).

What the sorting (i.e. starting at mode) does it makes the number of comparisons to make before stopping smaller, so it should be more efficient* then the "naive" algorithm using unsorted probabilities (ordinary inverse transform).

However my question is: does ignoring the sorting step have any other effect on the algorithm besides decreasing the performance of the algorithm (I guess not, but I'm looking for justification)? Moreover, how much (when) does it actually change in terms of performance?

_{* - I made some benchmark tests on this but they did not give any conclusive results.}

Answer 1

The algorithm recommended in the question is not guaranteed to be the most efficient in terms of comparisons (even ignoring the number of comparisons done during the sorting phase.) Consider the following algorithm:

1. (Initialization). Calculate the cumulative probabilities $p_k, \space k = 0,\dots, K$,
2. Sample $u$ from the Uniform distribution on $(0,1)$,
3. Find the smallest $k$ such that $p_k > u$ using binary search (https://en.wikipedia.org/wiki/Binary_search_algorithm).

This works because the cumulative probabilities are ordered, and it takes on average $\log_2 K-1$ comparisons to find the answer no matter what the distribution of the probabilities is. The worst-case performance is $\log_2 K + 1$ for the worst possible choice of $u$, regardless of distribution. The algorithm presented in the OP has a worst-case expected performance of $K/2$, achieved when all the $p_k$ are equal, but a best-case performance arbitrarily close to $1$ (when one of the $p_k$ has a probability very near $1$.) It can take as many as $K-1$ comparisons for any given $u$ no matter what the distribution is.

Note that for distributions with nonincreasing probabilities, such as the Geometric, the actual order of the elements of the vector being searched is the same for the two algorithms, and is the same as the order of the values themselves. The recommended algorithm becomes a linear search on the vector of cumulative probabilities, which makes for easy comparisons between the two. Interestingly enough, if the lower bound on the values is 0, the expected number of comparisons is equal to the expected value of the variate + 1! We can see that by observing that, in order to return a given value $x$, $x+1$ comparisons must have been made to reach that position in the list (the first position corresponding to $x=0$); and, given that the generator works, it must be that the probability of returning $x$ is equal to $p_x$. Therefore the probability of making $x+1$ comparisons is $p_x$, so the expected value of the number of comparisons is $\text{E}x+1$. Obvious alterations can be made for other lower bounds.

Clearly (in this case) if $\text{E}x+1 \geq \log_2K-1$ we would prefer the binary search algorithm. This comparison can obviously be extended to the general case by calculating the expected value of a variate on $[1, \dots, K]$ with probabilities $q_k$ corresponding to the sorted probabilities $p_j$ from the original distribution on whatever the range of $j$ was.

It would seem that the more heavily the distribution is concentrated in a small fraction of values, the better the recommended algorithm will be, but if you want something with guaranteed good performance both in terms of expected comparisons across distributions and worst-case comparisons, the binary search approach may be better. (Of course, if $K$ is small enough, a linear search on the CDF will outperform a binary search in runtime if not comparison count, as it avoids calculating and tracking the upper and lower bounds of the interval in which the answer is known to lie.)

Answer 2

Rather than expressing this algorithm in terms of search of a value satisfying an inequality, you can express it equivalently as a simple sum of indicators. Since you are interested in efficiency, it is also useful to vectorise the algorithm, allowing you to generate $m$ independent sample values. With this setup, the algorithm can be expressed as follows:

Vectorised sampling algorithm:

1. You have a probability vector $\boldsymbol{p} = (p_1,...,p_n)$.

2. Sort this by permutation $\pi$ to get $\boldsymbol{p}_* = (p_{\pi(1)},...,p_{\pi(n)})$. (If you don't want to sort then you can just take $\boldsymbol{p}_* =\boldsymbol{p}$.)

3. Calculate the cumulative permuted probabilities $(F_1,...,F_{n-1})$ given by: $$F_k \equiv \sum_{i=1}^{k} p_{\pi(i)}.$$

4. Simulate continuous uniform random variables $U_1,...,U_m \sim \text{IID U}(0,1)$.

5. Set $K_1,...,K_m$ as follows: $$\pi(K_i) \equiv 1 + \sum_{k=1}^{n-1} \mathbb{I}( U_i > F_k ).$$ Note: In the summation formula, once you get a single zero indicator, all later indicators are zero, so you can terminate the sum early, as soon as you get a zero indicator.

From this re-statement of the sampling algorithm, we can see that the benefit of sorting the probabilities into decreasing order is that the partial sums of the reordered probabilities start with the biggest values and so they increase rapidly early on. This makes it more likely that you will get to a zero indicator earlier, which allows you to terminate the sum, thereby saving some computational time.

However my question is: does ignoring the sorting step have any other effect on the algorithm besides decreasing the performance of the algorithm?

I can't see any other difference. Mathematically the algorithm works regardless of the permutation used (including no permutation). Thus, the only difference seems to be one of computational efficiency.

... how much (when) does it actually change in terms of performance?

Presumably you get an efficiency gain from sorting when the computational cost of applying and reversing the permutation $\pi$ is smaller than the computational gain from ending the summations earlier in the last step of the algorithm. Hence, you are more likely to get a larger computational gain when the number of categories $n$ is large and/or the number of samples $m$ is large. That is, the sorting step in the algorithm will tend to be more efficient when you are generating a larger number of samples from a categorical variable with a larger number of categories.

Answer 3

The ordering of the cumulative distribution is not needed, but its there automatically.

Let $(\phi_1,\ldots,\phi_n)$ be the categorical distribution and $p_i=\sum_{k=1}^i\phi_k$ the entries of the cumulative distribution. Now, since $p_1\le \ldots \le p_i\le p_{i+1}\le \ldots p_n$, the probability that a uniformly drawn $u\in (0,1)$ is in $(p_i,p_{i+1})$ is exactly $p_{i+1}-p_i=\phi_i$.

Even if the cumulative ordering is not ordered, that is you have $(p_{\pi(1)},\ldots,p_{\pi(n)})$, you would have to determine the $\pi(k)$ such that $p_{\pi(k)}=\max\{p_{\pi(i)}: p_{\pi(i)}\le u\}$ (max, because I sorted incresingly). Since the position of $p_{\pi(k)}$ in the unsorted cumulative distribution is $k$, you would chose $\phi_k$.