Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
1  
You should really ask support about this. If you do receive a reply, please post it! Also, thanks for the link to the alias method, which I didn't know. –  Szabolcs Mar 4 at 2:50
    
Regarding the performance: it's probably implemented in C, not in Mathematica, which might be the reason why it's faster than your Mathematica-based implementation. –  Szabolcs Mar 4 at 2:53
3  
You might find (if you've not already read it) wolfram.com/learningcenter/tutorialcollection/… interesting. Slightly dated but relevant, covers how MMA generally uses inverse CDF for such things until it deems that expensive, then switching to table-lookup or direct generation. As Szabolcs said, low-level internal code in C will always beat an interpreted interpretation... –  rasher Mar 4 at 3:05
    
@rasher +1 on your comment for "interpreted interpretation" ... clearly clear! –  belisarius Mar 4 at 3:17
    
@rasher Very interesting, it actually answered a related question I had today. Interesting enough that it's worth putting in a more prominent answer IMO. –  Szabolcs Mar 4 at 3:23

2 Answers 2

It inverts the CDF by binary search.

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.