[Numpy-discussion] random number generation in python compared to gsl

[Gabriel Gellner]

On Wed, Nov 05, 2008 at 03:19:09PM +0100, Matthieu Brucher wrote:
> > Not in this case: I always get the same sequence with seed=0
> > (different for both implementation, but the same each time I run it.)
> > I got around it by installing pygsl and taking random numbers from
> > there instead of from numpy.
> >
> > But I still find it strange to get two different sequences from two
> > implementation that claim to be the same algorithm...
> >
> > Giovanni
> 
> Hi,
> 
> I didn't check which MT was used, because there are several different.
> MT is based on one Mersenne prime, but for instance the Boost library
> wraps two generators with two different values. Perhaps the same
> occurs here.
> Other explanations include:
> - MT must use an array of starting values, perhaps the first is 0, but
> one is 0, 1, 2, 3, ... and the other uses 0, 0, 0, ...
This is the key, there is a lot of variation in the seeding algorithm that is
used in mersenne twister algorithms, usually they use some kind of linear
congruential algorithm to get started. I think that GSL uses the published
mersenne twister seed algorithm, what you can do is find out how gsl does this
and just set the full seed array in numpy.random.RandomState yourself. I have
done this for comparison with R's mersenne twister algorithm using my own MT
implementation and it works like a charm.

Gabriel