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
-------------------------------------------------------------