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

 
 

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

Submitted by:  Markus Bergholz <markuman>
Submitted on:  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

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( 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>
Project 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>
Project 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)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project 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.5