bugGNU Octave - Bugs: bug #42557, rand seed for twister algorithm...

 
 

bug #42557: rand seed for twister algorithm differs from Matlab

Submitter:  Markus Bergholz <markuman>
Submitted:  Sun 15 Jun 2014 08:28:31 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 19 Nov 2015 05:34:58 PM UTC, comment #7: 

Actually, scrap that last comment, 0 is seed 5489 in Matlab's latest RNG: http://www.mathworks.com/matlabcentral/answers/28216-how-to-duplicate-rand-twister-0-in-current-matlab#answer_36524

So only the normal distribution discrepancy remains.

Anonymous
Thu 19 Nov 2015 03:48:12 PM UTC, comment #6: 

The flipping of the last element of the state vector reference was supposed to link to bug #44046

Anonymous
Thu 19 Nov 2015 02:04:17 PM UTC, comment #5: 

This doesn't seem to hold for randn, tried with

Octave 4.0.0:

>> randn('twister', twister_seed(5489));
>> randn(1,5)

ans =

   1.362259   0.071959   0.203598  -0.096973   1.324298

Matlab R2015a:

>> rng(5489, 'twister');
>> randn(1,5)


ans =

    0.5377    1.8339   -2.2588    0.8622    0.3188

Also tried flipping the last state vector value as decribed (in both Matlab and Octave), randn still produces different streams.

Matlab also seems to use the Ziggurat method, so can't figure out where the difference is.

>> RandStream.getGlobalStream


ans =


mt19937ar random stream (current global stream)
             Seed: 5489
  NormalTransform: Ziggurat

Oddly enough, Matlab always seems to produce the same randn sequence no matter what the seed number is, i.e. rng(0, 'twister') and rng(5489, 'twister') generate the same 625 state vector when queried with s = rng

Anonymous
Sat 28 Jun 2014 05:15:20 AM UTC, comment #4: 


>> rand ('twister', [1 2 3]);
Error using rand
State must be a double scalar or the output of RAND('twister').

>> rng([1 2 3])
Error using rng (line 128)
First input must be a nonnegative integer seed less than 2^32, 'shuffle', 'default', or generator settings captured previously using S = RNG.


Markus Bergholz <markuman>
Thu 26 Jun 2014 03:41:21 PM UTC, comment #3: 

What does Matlab do if the seed is a vector, rather than a scalar?


rand ('twister', [1 2 3]);
rand



Rik <rik5>
Group administrator
Fri 20 Jun 2014 05:37:20 PM UTC, comment #2: 

and finaly


octave:8> rand('twister',twister_seed(5489))
octave:9> rand(5,1)
ans =

     0.814723686393179
     0.905791937075619
     0.126986816293506
     0.913375856139019
      0.63235924622541


seems to be the default seed from matlab...

Markus Bergholz <markuman>
Mon 16 Jun 2014 12:25:10 AM UTC, comment #1: 

I've marked the bug as confirmed and changed the category to "Matlab Compatibility".

Rik <rik5>
Group administrator
Sun 15 Jun 2014 08:28:31 AM UTC, original submission:  

Initial situation: 

Ruby, Python and Numpy using twister algorithm as default.
We can force this algorithm in Octave and Matlab too.
So, if you seed each software with the value 0, it should return the same random numbers.


Octave 3.8.1

octave:1> rand('twister',0)
octave:2> rand
ans =    0.844421851525048



Python 3.4.1

>>> import random
>>> random.seed(0)
>>> random.random()
0.8444218515250481



Python 3.4.1 Numpy 1.8.1

>>> import numpy
>>> numpy.random.seed(0)
>>> numpy.random.random()
0.5488135039273248




Ruby 2.1.2

irb(main):001:0> prng = Random.new(0)
=> #<Random:0x00000001d469e8>
irb(main):002:0> prng.rand()
=> 0.5488135039273248



Matlab 2013a

>> rand('twister',0)
>> rand

ans =

    0.548813503927325



So Octave is handling the seed input wrong. Afaiu the seed is processed in http://hg.octave.org/octave/file/03c2671493f9/liboctave/cruft/ranlib/getsd.f due use_old_generators. But maybe I'm wrong. However, This behavior can be fixed with another function.


function ret = twister_seed(SEED=0)

    ret = uint32(zeros(625,1));
    ret(1) = SEED;
    for N = 1:623
        ## initialize_generator
        # bit-xor (right shift by 30 bits)
        uint64(1812433253)*uint64(bitxor(ret(N),bitshift(ret(N),-30)))+N; # has to be uint64, otherwise in 4th iteration hit maximum of uint32!
        ret(N+1) = uint32(bitand(ans,uint64(intmax('uint32')))); # untempered numbers
    endfor
    ret(end) = 1;

endfunction


This function needs to be translated into c++ (that's past my c++ skills!) and need to be called in line 276 of libinterp/corefcn/rand.cc and returned to ColumnVector s (afaiu).
http://hg.octave.org/octave/file/03c2671493f9/libinterp/corefcn/rand.cc#l276



Octave 3.8.1

octave:1> rand('twister',twister_seed) # notice: default seed is 0 in the function
octave:2> rand
ans =    0.548813503927325
octave:3> rand
ans =    0.715189366372419
octave:4> rand
ans =    0.602763376071644
octave:5> rand('twister',twister_seed(42))
octave:6> rand
ans =    0.374540118847363
octave:7> rand
ans =    0.950714306409916
octave:8> rand
ans =    0.731993941811405




Matlab 2013a

>> rand('twister',0)
>> rand

ans =

   0.548813503927325

>> rand

ans =

   0.715189366372419

>> rand

ans =

   0.602763376071644

>> rand('twister',42)
>> rand

ans =

         0.374540118847362

>> rand

ans =

         0.950714306409916

>> rand

ans =

         0.731993941811405




Python 3.4.1 Numpy 1.8.1

>>> import numpy
>>> numpy.random.seed(0)
>>> numpy.random.random()
0.5488135039273248
>>> numpy.random.random()
0.7151893663724195
>>> numpy.random.random()
0.6027633760716439
>>> numpy.random.seed(42)
>>> numpy.random.random()
0.3745401188473625
>>> numpy.random.random()
0.9507143064099162
>>> numpy.random.random()
0.7319939418114051


Markus Bergholz <markuman>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #31549:  twister_seed.m added by markuman (456B - text/x-objcsrc)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by markuman (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-07-22 rik5 Dependencies- bugs #54342 is dependent
    2018-04-27 mtmiller Carbon-CopyRemoved 80942 -
    2018-04-27 mtmiller Release3.8.1 dev
    2014-06-16 rik5 Item GroupIncorrect Result Matlab Compatibility
        StatusNone Confirmed
        Summaryrand seed for twister algorithm rand seed for twister algorithm differs from Matlab
    2014-06-15 markuman Attached File- Added twister_seed.m, #31549

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code