7
$\begingroup$

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?

$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Mar 4, 2014 at 2:50
  • $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Mar 4, 2014 at 2:53
  • 4
    $\begingroup$ 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... $\endgroup$
    – ciao
    Mar 4, 2014 at 3:05
  • $\begingroup$ @rasher +1 on your comment for "interpreted interpretation" ... clearly clear! $\endgroup$ Mar 4, 2014 at 3:17
  • $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Mar 4, 2014 at 3:23

2 Answers 2

7
$\begingroup$

It inverts the CDF by binary search.

$\endgroup$
4
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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