[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
-------------------------------------------------------------
- rand("state"), Mike Miller, 2006/02/15
- Re: rand("state"), Bill Denney, 2006/02/15
- Re: rand("state"),
Mike Miller <=
- Re: rand("state"), Mike Miller, 2006/02/15
- Re: rand("state"), Paul Kienzle, 2006/02/15
- Re: rand("state"), Mike Miller, 2006/02/15
- Re: rand("state"), Mike Miller, 2006/02/15
- Re: rand("state"), Francesco Potorti`, 2006/02/16
- Re: rand("state"), Paul Kienzle, 2006/02/16
- Re: rand("state"), Mike Miller, 2006/02/16
- Re: rand("state"), Paul Kienzle, 2006/02/16
- Re: rand("state"), Francesco Potorti`, 2006/02/16
- Re: rand("state"), Mike Miller, 2006/02/16