python vs octave random generator

numpy, octave, python, random

Solution

Unfortunately the MT19937 generator in Octave does not allow you to initialise it using a single 32-bit integer as `np.random.RandomState(4)` does. If you use `rand("seed",4)` this actually switches to an earlier version of the PRNG used previously in Octave, which PRNG is not MT19937 at all, but rather the Fortran `RANDLIB`.

It is possible to get the same numbers in NumPy and Octave, but you have to hack around the random seed generation algorithm in Octave and write your own function to construct the state vector out of the initial 32-bit integer seed. I am not an Octave guru, but with several Internet searches on bit manipulation functions and integer classes in Octave/Matlab I was able to write the following crude script to implement the seeding:

function state = mtstate(seed)

state = uint32(zeros(625,1));

state(1) = uint32(seed);
for i=1:623,
   tmp = uint64(1812433253)*uint64(bitxor(state(i),bitshift(state(i),-30)))+i;
   state(i+1) = uint32(bitand(tmp,uint64(intmax('uint32'))));
end
state(625) = 1;

Use it like this:

octave:9> rand('state',mtstate(4));
octave:10> rand(1,5)
ans =

   0.96703   0.54723   0.97268   0.71482   0.69773

Just for comparison with NumPy:

>>> a = numpy.random.RandomState(4)
>>> a.rand(5)
array([ 0.96702984,  0.54723225,  0.97268436,  0.71481599,  0.69772882])

The numbers (or at least the first five of them) match.

Note that the default random number generator in Python, provided by the `random` module, is also MT19937, but it uses a different seeding algorithm, so `random.seed(4)` produces a completely different state vector and hence the PRN sequence is then different.

Problem

More specifically, numpy: ``` In [24]: a=np.random.RandomState(4) In [25]: a.rand() Out[25]: 0.9670298390136767 In [26]: a.get_state() Out[26]: ('MT19937', array([1248735455, ..., 1532921051], dtype=uint32), 2,0,0.0) ``` octave: ``` octave:17> rand('state',4) octave:18> rand() ans = 0.23605 octave:19> rand('seed',4) octave:20> rand() ans = 0.12852 ``` Octave claims to perform the same algorithm (Mersenne Twister with a period of 2^{19937-1}) Anybody know why the difference?

Original source