Alias Method, in depth: constant-time sampling from a discrete distribution, the Vose construction, exact integer tables and when a sum tree wins

The idea: fill n columns exactly

Start with n outcomes and probabilities p[i] that sum to 1. Multiply each by n, so the scaled values average exactly 1. Picture n columns, each of height 1, and n bars of heights n times p[i]. Some bars are shorter than 1 and some are taller. The claim at the heart of the method is that you can always cut the tall bars into pieces and use them to top up the short ones so that every column ends up exactly full, and that each column needs pieces from at most two outcomes: its own outcome at the bottom and one donor, the alias, on top.

If that holds, sampling is trivial. Pick a column uniformly at random, with probability 1 over n. Inside the column, flip a biased coin: with probability prob[col] return the column's own outcome, otherwise return alias[col]. The probability of any outcome j is the sum over columns of (1 over n) times the height of j's pieces in that column, and by construction those heights add up to n times p[j], so the result is exactly p[j].

Why can two pieces per column always be enough? Because the scaled values average 1, whenever some bar is shorter than 1 another bar must be taller than 1. Take one of each. The short bar fills the bottom of a column; the tall bar donates exactly the missing amount to fill it. That column is finished and never touched again. The tall bar is now shorter by the donated amount, and it goes back into the short pile or the tall pile depending on what is left. Each step finishes one column, so after n steps every column is full. This is an invariant argument: the remaining bars always average exactly 1 over the remaining columns, so the pairing can never get stuck.

The data structure

The data structure is two arrays of length n: prob, the height of the column's own outcome, and alias, the outcome that fills the rest. Nothing else is needed at sampling time, which is why the method is so fast and so cache-friendly: one draw touches one entry of each array.

Building the table in linear time

The Vose construction keeps two worklists, small and large, of indices whose scaled value is below 1 and at least 1. It repeatedly pairs one of each, exactly as in the argument above. Two details matter in floating point. The update of the large bar is written as (large + small) minus 1 rather than large minus (1 minus small), which loses less precision. And when one list empties, anything left in either list is set to probability 1: those values should be exactly 1 in exact arithmetic, and rounding is the only reason they are not.

import random
from collections import Counter

def build_alias(weights):
    n = len(weights)
    total = float(sum(weights))
    if n == 0 or total <= 0 or any(w < 0 for w in weights):
        raise ValueError("need non-negative weights with a positive sum")
    scaled = [w * n / total for w in weights]
    prob = [0.0] * n
    alias = list(range(n))
    small = [i for i, s in enumerate(scaled) if s < 1.0]
    large = [i for i, s in enumerate(scaled) if s >= 1.0]
    while small and large:
        s = small.pop()
        l = large.pop()
        prob[s] = scaled[s]                       # own share of column s
        alias[s] = l                              # l tops it up
        scaled[l] = (scaled[l] + scaled[s]) - 1.0 # what l has left
        (small if scaled[l] < 1.0 else large).append(l)
    for i in large + small:                       # leftovers are 1 up to rounding
        prob[i] = 1.0
    return prob, alias

def sample(prob, alias, rng=random):
    i = rng.randrange(len(prob))                  # uniform column
    return i if rng.random() < prob[i] else alias[i]

if __name__ == "__main__":
    w = [1, 2, 3, 4]
    prob, alias = build_alias(w)
    rng = random.Random(42)
    N = 1_000_000
    counts = Counter(sample(prob, alias, rng) for _ in range(N))
    for i, wi in enumerate(w):
        print(i, wi / sum(w), counts[i] / N)

Each index enters a worklist once at the start and re-enters at most once per pairing, and every pairing retires one index, so construction is O(n) time and O(n) memory. A zero weight is fine: its scaled value is 0, it goes to the small list, gets probability 0 and an alias, so its column always returns the alias.

Worked example: weights 1, 2, 3, 4

Trace the construction for weights 1, 2, 3 and 4, which are probabilities 0.1, 0.2, 0.3 and 0.4. Scaling by n = 4 gives 0.4, 0.8, 1.2 and 1.6. The small list starts as [0, 1] and the large list as [2, 3]; the code pops from the end of each.

Check the result by summing pieces. Outcome 0 appears only at the bottom of column 0: 0.25 times 0.4 = 0.1. Outcome 1: 0.25 times 0.8 = 0.2. Outcome 2 owns column 2 and tops up column 3: 0.25 times (1.0 + 0.2) = 0.3. Outcome 3 owns 0.8 of its own column and tops up columns 0 and 1 with 0.6 and 0.2: 0.25 times (0.8 + 0.6 + 0.2) = 0.4. Running the code above with a million draws reproduced all four probabilities to within 0.001. Note that the table is not unique; a different pop order gives a different but equally correct table.

Sampling millions at once

When you need millions of samples at once, as in a simulator or a graph-embedding pipeline, the per-draw Python overhead dominates. The sampling step vectorizes perfectly because each draw is independent and uses only two array lookups.

import numpy as np

def sample_many(prob, alias, size, rng):
    prob = np.asarray(prob)
    alias = np.asarray(alias)
    cols = rng.integers(0, len(prob), size=size)
    coins = rng.random(size)
    return np.where(coins < prob[cols], cols, alias[cols])

rng = np.random.default_rng(0)
prob, alias = build_alias([1, 2, 3, 4])
draws = sample_many(prob, alias, 10_000_000, rng)
print(np.bincount(draws) / draws.size)

A related trick uses a single uniform u in [0, 1): take the column as the integer part of u times n and the coin as the fractional part. It saves one random draw but the fractional part has fewer bits of precision as n grows, so prefer two independent draws unless the random number generator is the measured bottleneck.

Exact integer tables

Floating-point alias tables are accurate to roughly the precision of the stored probabilities, which is plenty for most simulation. When the weights are integers and exactness matters, for example in a test oracle, a lottery or a reproducible benchmark, build the table in integers. Scale every weight by n so that the column capacity is the total weight T, run the same pairing with integer arithmetic, and store integer thresholds. Sampling then draws a column with randrange(n) and a coin with randrange(T) and compares integers, so every outcome has probability exactly w[i] over T with no rounding anywhere. The cost is wider integers and a slightly more expensive coin.

Alias tables against the alternatives

The alias method is one of four standard ways to sample a discrete distribution. Which one is right depends on how many draws you take per build and whether the weights change.

The break-even is simple to reason about. With k draws per build, alias costs about a times n plus k and binary search about b times n plus k log n, where the build constant a is several times b because the alias build does more work per element than a prefix sum. Alias therefore wins once k exceeds roughly (a minus b) times n over log n, which in practice means a few times n draws or more. If weights change after every few draws, the O(n) rebuild dominates and a tree that updates one weight in O(log n) wins; prioritized replay buffers in reinforcement learning use a sum tree for exactly this reason. The tree structure is the one described in Fenwick Tree.

The same reasoning explains why language-model decoding does not build alias tables. The next-token distribution is different at every step, so each table would be used for one draw; a single pass over the probabilities is cheaper. Truncation methods such as top-p sampling already sort or partially sort the distribution, and sampling from the survivors is a short scan.

Where it is used

Good fits share one property: a fixed distribution sampled many times.

• Random walks on weighted graphs. Each node gets an alias table over its outgoing edge weights, built once, and every walk step is O(1). The node2vec paper describes precomputing alias tables for its biased transition probabilities for this reason. Memory is O(edges), which is the real limit on very large graphs.
• Weighted load balancing. Route requests to backends in proportion to capacity weights; rebuild only when the backend set or weights change.
• Simulation and synthetic data. Draw categories, event types or Zipf-like item popularity for millions of synthetic records from a fixed distribution.
• Negative sampling. Draw negatives from a fixed noise distribution, such as unigram counts raised to a power, many times per training step.

Failure modes

• Unnormalized or negative weights. Negative weights produce scaled values below zero and a table that silently returns garbage; validate inputs.
• Forgetting the leftovers. If the final loop that sets leftover entries to 1 is omitted, a column whose value drifted to 0.9999999 returns its alias, which may be an uninitialized default, with a tiny but real probability.
• Stale tables. The weights change but the table is not rebuilt, so the sampler follows yesterday's distribution. Version the table with the weights it came from.
• Rebuilding per draw. Code that builds a table and draws once has turned an O(n) scan into an O(n) build plus overhead.
• Testing by eye. Check every sampler with a chi-squared test against the target probabilities, with enough draws that the smallest expected count is in the hundreds; Probability in Algorithms covers the underlying bounds.

What to do next

1. Implement build_alias and sample from this page and run the empirical check on weights 1, 2, 3, 4.
2. Add a chi-squared test for a random 1,000-outcome distribution with ten million draws.
3. Profile your current sampler: count draws per rebuild and time both builds to find the break-even k for your n.
4. If weights change often, switch to a sum tree; if they are static, switch to an alias table.
5. For integer weights where exactness matters, implement the integer variant and test it with exact expected counts.
6. For a streaming source where you cannot hold all weights, use reservoir sampling instead.