numpy performance and random numbers

[Robert Kern]

On 2009-12-19 09:14 AM, Carl Johan Rehn wrote:
> On Dec 19, 2:49 pm, sturlamolden<sturlamol... at yahoo.no>  wrote:
>> On 19 Des, 11:05, Carl Johan Rehn<car... at gmail.com>  wrote:
>>
>>> I plan to port a Monte Carlo engine from Matlab to Python. However,
>>> when I timed randn(N1, N2) in Python and compared it with Matlab's
>>> randn, Matlab came out as a clear winner with a speedup of 3-4 times.
>>> This was truly disappointing. I ran tthis test on a Win32 machine and
>>> without the Atlas library.
>>
>> This is due to the algorithm. Matlab is using Marsaglia's ziggurat
>> method. Is is the fastest there is for normal and gamma random
>> variates. NumPy uses the Mersenne Twister to produce uniform random
>> deviates, and then applies trancendental functions to transform to the
>> normal distribution. Marsaglia's C code for ziggurat is freely
>> available, so you can compile it yourself and call from ctypes, Cython
>> or f2py.
>>
>> The PRNG does not use BLAS/ATLAS.
>
> Thank you, this was very informative. I know about the Mersenne
> Twister, but had no idea about Marsaglia's ziggurat method. I will
> definitely try f2py or cython.

It's worth noting that the ziggurat method can be implemented incorrectly, and 
requires some testing before I will accept such in numpy. That's why we don't 
use the ziggurat method currently. C.f.

   http://www.doornik.com/research/ziggurat.pdf

Requests for the ziggurat method come up occasionally on the numpy list, but so 
far no one has shown code or test results. Charles Harris and Bruce Carneal seem 
to have gotten closest. You can search numpy-discussion's archive for their 
email addresses and previous threads. Additionally, you should ask future numpy 
questions on the numpy-discussion mailing list.

   http://www.scipy.org/Mailing_Lists

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco