numpy performance and random numbers

[sturlamolden]

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.