numpy performance and random numbers

[Carl Johan Rehn]

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.

Well, I guess I knew that random numbers were not handled by BLAS/
ATLAS, but wasn't 100% sure.

Carl