For some reason alias sampling is not well known, even though it strictly dominates using a search tree when generating randomness (it's faster, it uses less memory, and IMO it's easier to code). It even works better when generalized to the case of preparing superpositions on quantum computers ( https://algassert.com/post/1805 ).

Alias sampling is a really beautiful algorithm, but implementing it correctly and efficiently isn't as easy as you seem to imply.

I would also very much like to know what a generalization of this algorithm to continuous distribution would look like.

You make a list of below-average bars and a list of above-average bars. Iteratively pop the top of both, donate from the larger to the smaller until the smaller is average, then put the larger back into the appropriate list based on its new value. Stop when there's no too-large bars left. Compared to writing a red black tree it's trivial.

A correct implementation wrt probabilities being exact given a frequency table (the "average" can be a very large rational number), as well as efficient when constructing the data structure to be sampled from during sampling phase is IMO not as trivial as you seem to imply.

I have asked people to implement this algorithm during engineering interviews on a regular basis for the better part of five years, and I can promise you, strictly none of the candidate who were subjected to it would describe their experience as "I was asked a trivial question".

But then everyone has their own bar for the word "trivial". I know for example that I have my own bar for the word "arrogance".

Oh wow, how do you prompt the interviewees? Do you tell them about the average value threshold or the guarantee you only need 2N buckets? If not, I’d guess this is probably way too hard of an interview question, even though a decent implementation really isn’t that difficult. Is this like a bonus difficult question for someone who aces all your other questions?

The implementation isn’t hard once you’ve read the Vose paper or the Keith Schwarz link I posted. There’s only one tricky part about handling rounding that is pretty easy to fix.

I would never expect anyone who hasn’t heard of the idea to be able to implement it on the spot without being more or less given a recipe for the math and how to use it. If you ever find someone who can, hire them immediately!

This is one of those problems that only feels somewhat easy in retrospect, but almost nobody could derive any of it on their own. The alias method is extremely clever, and not at all trivial.

> Oh wow, how do you prompt the interviewees?

Back when I was still doing this, I used to start the interview like this:

        - given: a non-uniform, finite, frequency table with actual observations of some real-world phenomenon (and if that's too generic, just tell the candidate you've thrown a dice a million time an recorded the outcomes).
        - given: a uniform pseudo random number generator
        - please design a pseudo random number generator that will, if sampled for long enough, produce the same distribution (non generic version: RNG that statistically behave like the dice).

Unfortunately, a non-negligible proportion of candidates fail to produce any solution at all :(

Decent candidates quickly produce the standard: "build the CDF once and then at sample time, draw a uniform number in the correct range and walk the CDF until we're above that random number, return the bucket where that occurred"

I then ask for time and space complexity for the sampling phase.

Slightly better than average candidate notice that the CDF is naturally sorted and can be binary searched.

I then ask those if they can improve what we have to a O(1) solution (independent of the size of the histogram).

The top 20% candidates usually manage to put their intuition at work: from the fact that O(1) almost always implies a simple table lookup, they come up with the: "build a very large array, where the index of each bucket B of the histogram is replicated histo[B] times and then sample uniformly from it".

We then discuss when this algorithm is/isn't practical.

The candidates that get to this point within ~15mn can then be gently walked towards redesigning the alias method using some hints.

The starting hint goes like this: what if you were to completely ignore the fact that the histogram is not uniform, and just stubbornly sample from it, can you quantify how much of a mistake you would make ?

The second hint is: for some buckets, picking them will be an overshoot, and for some we will undershoot. Explain when that is.

The third hint: say we picked a bucket that's chosen too often, can we decide to only keep it some of the time (as in: with some carefully chosen probability)

The fourth hint: and if we decide not keep it, why waste the dice throw, why not pick another, carefully chosen bucket

etc ...

The nice property of this question is: it's very gradual, starting with a problem most people with a CS degree should be able to solve, then going to a solution that's less standard and requires more than base regurgitation of CS lessons, to finally get to what I consider to be (at the risk of repeating myself) the non-trivial state of the art solution to the problem.

The point where the candidate start to struggle gives you a nice reading on their strength.

And there's the rare few new grads that would vaguely recall having seen the alias method and could rebuild it cleanly on the white board from scratch and without help. These rare occurrences were very enjoyable.

I agree that the discretizing step is hard to get right. I agree that you don't have enough time in an engineering interview to get Vose's algorithm right (nevermind think of it in the first place). I was definitely imagining that you had a day to implement it and test it, not an hour and a whiteboard.

The closest continuous distribution analogy would be inverse transform sampling - analytically determine the CDF and invert it. Then you have an O(1) function to map a uniform sample into your continuous PDF. Unfortunately, it’s not always possible, lots of continuous functions can’t be integrated & inverted.

There's actually an interesting experiment to run, which is as follows:

       - take some non-uniform a continuous distribution with finite support.
       - discretize it into a N-bucket histogram
       - apply the alias algorithm
       - display the result
       - increase N
       - rinse and repeat

The resulting animation is likely going to be interesting.

I suspect it won't be stable at all and might veer into the fractal at some point.

Should be easy to try, maybe I’ll do it if I understand correctly what you’re thinking about...

It sounds like you’re thinking of a feedback loop, where I pick, say 1000 uniform samples, then run them through the alias map, and then run those results back through the alias map, etc., while each time increasing the number of buckets for the alias map?

I would assume the output distribution range is either [0,1] or we re-normalize the samples every time through the alias map?

Is that what you’re thinking of, or are you suggesting to visualize the map directly, and not a bunch of samples?

For a feedback loop on samples, I guess I might expect the samples to eventually converge with all samples at the function’s maximum. But if you mean to track and visualize all the little intervals in the alias map through itself several times, yeah, that might be a bit wild to see.

Thanks for the links. The last one is long, but made it easy to understand how and why alias method works (and it's awesome for providing mathematical proofs along with algorithms).

That said, I'm having trouble understanding the solution you described. I can understand how the alias method can help me efficiently sample an arbitrary discrete probability distribution, but I can't see how adding jitter helps with the problem from TFA - namely, to generate samples that conform to a given distribution long-term, but avoid hitting "annoying" low-probability outcomes like strings of consecutive lower-probability options. Could you please clarify / expand on your explanation?

Sure, my thinking is that N stratified uniform samples plugged into the alias map is the part that prevents the annoying low-probability sequences. This is very similar to the simpler and more common technique of picking from a shuffled list of outcomes, but plugging it into the alias map allows the outcomes to have arbitrary probabilities not quantized by the number of choices.

The shuffle and jitter are the things that control randomness. Shuffle controls the order of events, and jitter controls which events are chosen. Using both you can have different outcomes each time through the N samples.

Regular sampling would give you the same set of events, where jittering will change which events are chosen depending on your parameter choices. But! Jitter will only change the number of each event chosen by a limited amount that you control.

The choices for the blending knob and the number N of uniform samples, along with your goal distribution, let you control how often and by how much the events change. For example, by picking an N that gives you uniform strata that are the same size as the smallest item in the goal distribution, you can guarantee that you get at least one choice from each item in the goal distribution. I think this would also ensure that adding jitter changes the number of choices for each goal bucket by at most 1. It might be by at most 2, due to the alias map’s splitting of large buckets, and this cap might need strata that are half the size of the smallest bucket. I’m not sure of the exact details without trying it, but I am pretty sure this gives you control over the amount of variation; the number of events you can possibly have above or below the expected value.

That’s the idea though - choosing N and using the knob allows you to decide how much variation you want on top of the guarantee that you never get sequences that feel wholly improbable.

Thanks for elaborating!

I did some related things like implementing stratified sampling for RNG used in Monte Carlo simulations, so I understand some of the concepts here, but I still don't have a good intuition for how jitter and shuffling affects the outcome set, or whether it could replicate manually preventing repetitions in results from happening. However, you provided enough information that I think I can figure it out now.

Tell me if I'm understanding this correctly:

- You generate a sequence of uniform random numbers using stratified sampling. This increases variation by ensuring the [0, 1) space is sampled more-less evenly, while not breaking the uniformness of resulting distribution.

- With strata size set to the lowest probability in the desired outcome set, we can ensure even the lowest probability outcome will happen n times (with n = samples per stratum). This allows the low probability events to be observed like a naive expectation would suggest, while still not breaking the overall distribution.

- Alias map is used to translate the stratified uniform sequence to desired distribution (is this the discrete equivalent of rejection sampling?).

- Shuffling is applied to the stratified uniform sequence, and ensures events don't happen in order of decreasing probabilities (modulo shuffling introduced by alias map breaking probabilities into pieces).

- Jitter lets us add additional slight unpredictability, while still not breaking the desired distribution (because jitter is uniformly distributed, so it averages out).

Do I understand this correctly?

Yes, spot on, I think that describes it exactly.

And yes, the alias map serves a very similar function to rejection sampling of a discrete distribution, except that you don’t have to waste any samples. And the alias map serves as a deterministic mapping function; the same input will always yield the same output.

There’s another method for discrete distribution sampling that is closer to alias mapping, and better preserves stratification than the alias method but costs O(log n) per sample: compute the discrete CDF, and then invert it by picking a uniform sample and using that to binary search on the (sorted by definition) CDF to recover the choice, which will be then distributed according to your original discrete PDF array. For sample set size that is known in advance, you can avoid the binary search, and do this in O(1) time per sample.

The method in the article feels similar-ish to the binary search method, maybe it’s equivalent, but I’m not sure.

Nice, seems like a great solution especially in the use case where this service is needed by a game server and good performance is a requirement.