Alias method

In computing, the alias method is a family of efficient algorithms for sampling from a discrete probability distribution, due to A. J. Walker.^[1][2] That is, it returns integer values 1 ≤ i ≤ n according to some arbitrary probability distribution p_i. The algorithms typically use O(n log n) or O(n) preprocessing time, after which random values can be drawn from the distribution in O(1) time.^[3]

Operation

Internally, the algorithm consults two tables, a probability table U_i and an alias table K_i (for 1 ≤ i ≤ n). To generate a random outcome, a fair die is rolled to determine an index into the two tables. Based on the probability stored at that index, a biased coin is then flipped, and the outcome of the flip is used to choose between a result of i and K_i.^[4]

More concretely, the algorithm operates as follows:

1. Generate a uniform random variate 0 ≤ x < 1.
2. Let i = ⌊nx⌋ + 1 and y = nx + 1 − i. (This makes i uniformly distributed on {1, 2, …, n} and y uniformly distributed on [0, 1).)
3. If y < U_i, return i. This is the biased coin flip.
4. Otherwise, return K_i.

An alternative formulation of the probability table, proposed by Marsaglia et. al.^[5] as the “square histogram” method, uses the condition x < V_i in the third step (where V_i = (U_i + i − 1)/n) instead of computing y.

Table generation

The distribution may be padded with additional probabilities p_i = 0 to increase n to a convenient value, such as a power of two.

To generate the table, first initialize U_i = np_i. While doing this, divide the table entries into three categories:

• The “overfull” group, where U_i > 1,
• The “underfull” group, where U_i < 1 and K_i has not been initialized, and
• The “exactly full” group, where U_i = 1 or K_i has been initialized.

If U_i = 1, the corresponding value K_i will never be consulted and is unimportant, but a value of K_i = i is sensible.

As long as not all table entries are exactly full, repeat the following steps:

1. Arbitrarily choose an overfull entry U_i > 1 and an underfull entry U_j < 1. (If one of these exists, the other must, as well.)
2. Allocate the unused space in entry j to outcome i, by setting K_j = i.
3. Remove the allocated space from entry i by changing U_i = U_i − (1 − U_j) = U_i + U_j − 1.
4. Entry j is now exactly full.
5. Assign entry i to the appropriate category based on the new value of U_i.

Each iteration moves at least one entry to the “exactly full” category (and the last moves two), so the procedure is guaranteed to terminate after at most n−1 iterations. Each iteration can be done in O(1) time, so the table can be set up in O(n) time.

Vose^[3]:974 points out that floating-point rounding errors may cause the guarantee referred to in step 1 to be violated. If one category empties before the other, the remaining entries may have U_i set to 1 with negligible error.

As the lookup procedure is slightly faster if y < U_i (because K_i does not need to be consulted), one goal during table generation is to maximize the sum of the U_i. Doing this optimally turns out to be NP hard,^[5]:6 but a “Robin Hood” heuristic comes reasonably close: rob from the richest and give to the poorest. That is, at each step choose the largest U_i and the smallest U_j. Because this requires sorting the U_i, it requires O(n log n) time.

Efficiency

Although the alias method is very efficient if generating a uniform deviate is itself fast, there are cases where it is far from optimal in terms of random bit usage. This is because it uses a full-precision random variate x each time, even when only a few random bits are needed.

One case arises when the probabilities are particularly well balanced, so many U_i = 1 and K_i is not needed. Generating y is a waste of time. For example if p₁ = p₂ =  ¹⁄₂, then a 32-bit random variate x could be used to make 32 choices, but the alias method will only generate one.

Another case arises when the probabilities are strongly unbalanced, so many U_i ≈ 0. For example if p₁ = 0.999 and p₂ = 0.001, then the great majority of the time, only a few random bits are required to determine that case 1 applies.

In such cases, the table method described by Marsaglia et al.^[5]:1–4 is more efficient.

Literature

Knuth, Art of Computer Programming, Vol 2: Seminumerical Algorithms: Sect. 3.4.1.

Implementations

• http://www.keithschwarz.com/darts-dice-coins/ Keith Schwarz: Detailed explanation, numerically stable version of Vose’s algorithm, and link to Java implementation
• http://apps.jcns.fz-juelich.de/ransampl Joachim Wuttke: Implementation as a small C library.
• http://oroboro.com/non-uniform-random-numbers Rafael Baptista’s Implementation in C++

References