help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: rand("state")


From: Mike Miller
Subject: Re: rand("state")
Date: Wed, 15 Feb 2006 13:41:12 -0600 (CST)

On Wed, 15 Feb 2006, Bill Denney wrote:

On Wed, 15 Feb 2006, Mike Miller wrote:

On my system, rand("state"), using Octave-Forge rand.oct, returns a row vector of 625 integers. (By the way, the "help rand" doc says that it is a column vector.) The interesting thing is that the integers seem to be uniformly distributed with a max of about 2^32, but the 625th integer seems to be different from the others with a mean somewhere around 300 or so. What is rand("state")(625) doing?

If I remember correctly, the way that the random number generator works is with a Mersenne Twister algorithm. This algorithm uses what is called an incomplete matrix calculation and one of the rows of the matrix is shorter than others. I think that's why the last number is smaller.

OK. I knew it was Mersenne Twister, and I looked at some MT web pages but I didn't find what I needed. You suggestions seem reasonable.


stored_seeds = ceil( rand( 100, 625 ) * 2^32 );
rand ("state", stored_seeds( rep, : ) )

Where "rep" is a scalar holding the replicate number.

This is a question for someone else, but my suggestion would be to try out seeding with one seed and see if the numbers come out the same every time. My guess would be that the seeding algorithm would just chop off any extra bits that are unnecessary in the last number.

Using the method above, use of the same seed produced the same sequence, but the seed is altered somehow. I'm not sure that it is a problem for me, but this is what I'm seeing:

octave:165> stored_seeds = ceil( rand( 100, 625 ) * 2^32 );
octave:166> rep = 20; rand ("state", stored_seeds( rep, : ) );
octave:167> rand("state")(625)
ans = 1
octave:168> rand("state")(1)
ans = 2147483648
octave:169> stored_seeds (20,[1,625])
ans =

  1948246809  4267933407

The final (625th) integer in the seed always equals 1 even if the seed I request has a different value in that position.

What I am finding is that if the 625th integer is outside the range from 1 to 624, it is reset to 1 and all the other integers are changed, but if I set that final 625th integer to be within the range from 1 to 624, everything works correctly:

octave:171> rep = 20; rand ("state", [stored_seeds( rep, 1:624 ), 27] );
octave:172> rand("state")(1)
ans = 1948246809
octave:173> rand("state")(625)
ans = 27

So I think I'll just do this as it seems to work out fine:

octave:183> stored_seeds = ceil( [ rand( 100, 624 ) * 2^32 , rand( 100, 1 ) * 
624 ] );
octave:184> rep = 20; rand ("state", stored_seeds( rep, : ) );
octave:185> sum ( rand("state") == stored_seeds( 20, :) )
ans = 625


To summarize: I think this produces a good set of 100 random seeds for the MT generator:

stored_seeds = ceil( [ rand( 100, 624 ) * 2^32 , rand( 100, 1 ) * 624 ] );

That is a matrix of integers of size 100 x 625.

Mike



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]