Algorithm in RandomChoice[{w1,w2,...}->{e1,e2,...},n]?

Question

I wonder if someone happens to know the actual algorithm implemented in the particular overload of RandomChoice that handles explicitly presented discrete distributions, namely

RandomChoice[{w1, w2, ...}->{e1, e1, ...}]
    gives a pseudorandom choice weighted by the w1.

It's blazing fast, as illustrated in my test notebook here: https://www.dropbox.com/s/rwm80fut60v926b/FastNonUniformPseudoRandoms002.cdf

(btw, such things are really useful in, say, particle filters).

I know of at least four algorithms for solving this problem. A greedy-space array of outcomes has O(1) time but O(S) space where S is the sum of the frequency counts (and only works when the weights are integers or can be scaled up to integers). Two algorithms explicitly invert the CDF, either by linear search -- O(N) time, O(N) space, or binary search -- O(log N) time, O(N) space. The best I know is Walker's Alias Method http://en.wikipedia.org/wiki/Alias_method

which is O(N) space, O(1) time. I implemented this in my notebook and it beats the pants off the others (as expected) but is still 30 times slower than the built-in (though still just O(1)).

I wonder if the built-in is just an optimized Walker's or is it something better?

Answer 1

It inverts the CDF by binary search.

Answer 2

You might find (if you've not already read it) Random Number Generation from the Wolfram tutorial collection interesting. Slightly dated but relevant, covers how Mathematica generally uses inverse CDF for such things until it deems that expensive, then switching to table-lookup or direct generation. Probably worth a ping to support, they seem pretty open about most internals, and I would not be surprised if current versions use something like Marsaglia's fast methods. As Szabolcs said, low-level internal code in C will always beat an interpreted interpretation.