random seed for statistics toolbox
Edmund Brekke
Previously I have used the syntax
rand('state',randState);
randn('state',randState);
but now Matlab tells me that this is no longer the right way of doing it (by the way, that also gave some very strange results including a consistent positive bias for Gaussian numbers).
I therefore try to use the new syntax
stream=RandStream('mrg32k3a','Seed',randState);
to for example draw Rayleigh distributed random numbers:
raylrnd(stream,sigma)
This gives the following error message:
Warning: Size vector should be a row vector with integer elements.
??? Error using ==> RandStream.RandStream>throwUndefinedError at 756
Undefined function or method 'lt' for input arguments of type 'RandStream'.
Error in ==> RandStream.RandStream>RandStream.lt at 718
function a = lt(varargin), throwUndefinedError; end
Error in ==> raylrnd at 32
b(b < 0) = NaN;
Can anybody please help me understand what is going on? There must be some way of controlling the random seed for the statistics toolbox, but I cannot find any instructions on this in the help files.
Wayne King
Hi Edmund, You don't want to pass the argument stream to raylrnd(), you want to create the random number stream and then call raylrnd(). For example, if you want to produce consistent results from two calls to raylrnd(), you can do the following:
S = RandStream('mt19937ar');
RandStream.setDefaultStream(S);
R = raylrnd(2,10,1); % a 10X1 vector parameter =2
S = RandStream('mt19937ar');
RandStream.setDefaultStream(S);
R1 = raylrnd(2,10,1);
% compare R and R1
Hope that helps,
Wayne
Peter Perkins
> I have a problem with controlling the random seed.
> Previously I have used the syntax
>
> rand('state',randState);
> randn('state',randState);
>
> but now Matlab tells me that this is no longer the right way of doing it
It does, but it's also true that ALL OF THE OLD SYNTAX CONTINUES TO WORK AS IT ALWAYS HAS. So no code is broken. Can't stress that enough.
> (by the way, that also gave some very strange results including a
> consistent positive bias for Gaussian numbers).
You'd have to demonstrate that. I don't know of any such problem unless perhaps you are passing in something very strange as the second arg. By the way, unless you mistyped, passing the same "randState", in the sense of a full state vector, into both rand and randn would be an error. Perhaps you are really passing in a scalar integer.
> I therefore try to use the new syntax
> stream=RandStream('mrg32k3a','Seed',randState);
>
> to for example draw Rayleigh distributed random numbers:
>
> raylrnd(stream,sigma)
>
> This gives the following error message:
[snip]
As Wayne points out, RAYLRND is not prepared to accept a random number stream as an input in the same way that (the newest versions of) RAND and RANDN are. Let's figure out what it is you're trying to do. What follows seems long, but that's partially because I'm not sure what you need.
The old syntaxes above did two things:
1) they made sure a particular algorithm was "active" for the two separate random number generators behind RAND and RANDN.. The "subtract with borrow" algorithm in the RAND case, "congruential+shiftregister+ziggurat" in the RANDN case.
2) they set the states for the two generators.
I can't tell if you mean for "randState" to be an integer, in which case you're "initializing" the generators, or if it's a full state vector, in which case you're "restoring a previous state". If the latter, I'll assume you really meant
> rand('state',randState);
> randn('state',randnState); % randnState, not randState
i.e., two different state vectors.
Now the new syntax: First of all, RAND and RANDN no longer have separate generators behind them, so there's only one thing to control. Simpler. That thing is known as the "default stream", where default is in the sense of "the stream that RAND or RANDN use unless you say otherwise." All of the functions in the Statistics Toolbox, like RAYLRND, eventually call RAND or RANDN, so there's still only one thing to control for all of those too.
When you start up MATLAB, the default stream uses the Mersenne Twister algorithm (that's been true for a while now, since R2007b I think). That's a good choice for most purposes, so you probably don't need to create a new stream (unless you want multiple streams, which I'm guessing you don't). So it's just a matter of controlling the default stream that's already there. The help for RAND, RANDN, and RANDSTREAM have examples of what you need to do.
For example, if you want to be able to repeat a particular calculation, you want to save and restore the default stream's state. Here's the appropriate example in the help for RAND:
>> help RandStream
[snip]
Save and restore the default stream's state to reproduce the output of
RAND:
defaultStream = RandStream.getDefaultStream;
savedState = defaultStream.State;
u1 = rand(1,5)
defaultStream.State = savedState;
u2 = rand(1,5) % contains exactly the same values as u1
This is very similar to what you would have done with the old syntax using a full state vector as the second input argument. The main difference is that there's a line of code that gets the default stream explicitly. More typing, but less ambiguous about what's happening.
You may just want to initialize the random number generator, i.e., what you were doing by with the old syntax using an integer as the second input argument. That's somewhat simpler.
>> help RandStream
[snip]
Reset the generator for the default stream that underlies RAND, RANDI,
and RANDN back to the beginning to reproduce previous results:
reset(RandStream.getDefaultStream);
Again, more typing, but less ambiguous about what's happening. It's also possible to pass an integer into RESET as a seed, i.e.,
>> help RESET
[snip]
RESET(S,SEED) resets the generator for the random stream S to the initial
internal state corresponding to the seed SEED. Resetting a stream's seed
can invalidate independence with other streams.
Resetting a stream should be used primarily for reproducing results.
reset(RandStream.getDefaultStream,1234);
but as the help says, you may need to be careful when doing that if you are using multiple streams.
Hope this helps. I wrote a couple of entries a while back for Loren Shure's blog that talk about all this and also cover cases that are way more complicated. They may be of interest.
<http://blogs.mathworks.com/loren/2008/11/05/new-ways-with-random-numbers-part-i/>
<http://blogs.mathworks.com/loren/2008/11/13/new-ways-with-random-numbers-part-ii/>
Peter Perkins
> Hi Edmund, You don't want to pass the argument stream to raylrnd(), you
> want to create the random number stream and then call raylrnd(). For
> example, if you want to produce consistent results from two calls to
> raylrnd(), you can do the following:
>
> S = RandStream('mt19937ar');
> RandStream.setDefaultStream(S);
> R = raylrnd(2,10,1); % a 10X1 vector parameter =2
> S = RandStream('mt19937ar');
> RandStream.setDefaultStream(S);
> R1 = raylrnd(2,10,1); % compare R and R1
Wayne, nothing you have there is wrong per se, but it is probably more complicated than it needs to be:
1) MATLAB already starts up using the Mersenne Twister algorithm ('mt19937ar', to be precise) as the "default stream". So unless you want something else, there's no need to create a new stream and make it the default.
2) If you want to "start over", it's usually simpler to reset an existing stream rather than creating a new one.
So here's how I would have done the above:
S = RandStream.getDefaultStream();
reset(S); % leave this out if you're in a new MATLAB session
R = raylrnd(2,10,1); % a 10X1 vector parameter =2
reset(S);
R1 = raylrnd(2,10,1); % compare R and R1
More generally, you can "reseed", to get results that are different than what MATLAB gave you at startup:
S = RandStream.getDefaultStream();
reset(S,j);
R = raylrnd(2,10,1); % a 10X1 vector parameter =2
reset(S,j);
R1 = raylrnd(2,10,1); % compare R and R1
Another possibility is
S = RandStream.getDefaultStream();
savedState = S.State;
R = raylrnd(2,10,1); % a 10X1 vector parameter =2
S.State = savedState;
R1 = raylrnd(2,10,1); % compare R and R1
As the name "reset" implies, the first two are "reinitialization", a quick way to start over from the beginning. The third is more for "repeating specific calculations", it allows you to reproduce values from any point in the stream of random values (as long as you saved the state at the right point). Which one is "right" depends on what you are doing.
Hope this helps.