Prev: How to run Matlab 6.5 in Windows XP x64 Edition
Next: nested simulation or time-generating "for"
From: Peter Perkins on 16 Mar 2010 17:26
On 3/16/2010 2:40 AM, Kevin wrote:

> (2) This point is about the documentation (I am using R2009b) on whether
> rand and randn are using the same random number stream or not. Now after
> thinking about it more and reading your reply, the answer is yes (if we
> are using the new RandStream object. this is true for both default
> stream and user-created stream) but the answer is no (if we are using
> legacy random stream).

We're looking at what improvements the doc needs. Explaining the difference between "legacy mode" and the current behavior is definitely one of those areas. Thanks for pointing out these specific places we need to focus on.


> (3) This is about your last posting on "Save and restore state to
> reproduce results". Even if I save the state, I still cannot guarantee
> that I will get the same random numbers. For example, after I save the
> state, someone's code changes the default stream to use different
> algorithm, eg.
>
>>> dfltStream = RandStream.getDefaultStream;
>>> savedState = dfltStream.State;
>>> h = RandStream('mt19937ar', 'Seed', 101);
>>> RandStream.setDefaultStream(h)
>
> Things are not that simple anymore since Matlab allows user to set the
> default stream to something else. This is a gotcha. Correct me if I am
> wrong on this.

You are right, but that has ALWAYS been true, even when you were typing things like rand('state',0). And in that older world, there was no good way to even tell that that had happened. Now at least you can get the (current) default stream and easily tell if it's been changed out from under you.

Actually, if you run your example in a fresh MATLAB session, you _will_ get what you want, but not for the right reasons. If however, someone had typed 'mrg32k3a', then your saved state would not be relevant for the (now) current default stream. But again, that's ALWAYS been true.

I can't tell you what to do about this, because I don't know the context in which you are writing this code. Code that wants to be a "good citizen" should _never_ replace the default stream, any more than code buried in some function, and written for MATLAB prior to R2009b, should have done things like rand('state',0), or for that matter, "close all force". It has global side effects. There aren't many things in MATLAB that have global side effects but there are some. Messing with the default stream (other than drawing values from it) is something that code in functions should leave to the "end user" as much as possible, so that they know and control what is going on. I touch on this in the second of these two blog entries:

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




Just to elaborate on what would happen in your example: if your code did this:

dfltStream = RandStream.getDefaultStream;
savedState = dfltStream.State;

and then some code that you called did this:

RandStream.setDefaultStream(RandStream('mrg32k3a','Seed',101))

and then your code did this:

dfltStream.State = savedState;

or even

dfltStream = RandStream.getDefaultStream;
dfltStream.State = savedState;

you'd get an error, because the internal states for those two generators are a completely different format.

From: Kevin on 17 Mar 2010 01:10
Hi Peter,

Thanks for you reply to my question (4). Now I understand what is going on. Thanks for sharing the info here. Is it possible for you to mention this in the Matlab documentation, something like this: The following two methods will generate identical random streams except for seed = 0,

Method 1:
---------------
h = RandStream('mt19937ar', 'Seed', seed);
rand(h, 1, 3)

Method 2:
---------------
rand('twister', 0);
rand(1,3)

Now I think about this more. I get a bit confused. Do both RandStream('mt19937ar') and rand('twister', seed) use the same Mersenne twister algorithm? Or it is the mapping of 0 being different in the old and new ways?

This is a concern to me because I would like to replace the old rand and randn with RandStream and at the same time I want to make sure that I get the same random numbers.

I just happen to choose 0 as my seed in my example. I can add a checking in my code to avoid using seed 0.



Peter Perkins <Peter.Perkins(a)MathRemoveThisWorks.com> wrote in message <hnorhg$5a8$1(a)fred.mathworks.com>...
> On 3/16/2010 2:40 AM, Kevin wrote:
>
> > (4) I have posted this question in a new thread but I did not get any
> > reply. So I will post it here. I got the impression from the
> > documentation ("Legacy Mode") that I can use RandStream to generate the
> > same random numbers as generated by the legacy code. But I did the
> > following experiment and found out that I am wrong.
> >
> >>> h = RandStream('mt19937ar', 'Seed', 0); rand(h, 1, 3)
> >
> > ans =
> >
> > 0.8147 0.9058 0.1270
> >
> >>> rand('twister', 0); rand(1,3)
> >
> > ans =
> >
> > 0.5488 0.7152 0.6028
> >
> > So am I having the wrong concept thinking that I can use RandStream to
> > generate the same random numbers as generated by the legacy code.
>
> You've discovered the one inconsistency that I know of between old and new (other than rand and randn working separately). This is a bit arcane, and much/all of what follows is not something that anyone would really need to worry about, unless they happened to type exactly what you typed and said, huh?"
>
> You will note that none of the old help recommended using 0 as in rand('twister',0), it used 5489, which is the default seed used by the designers of the Mersenne Twister algorithm. It was not clear if 0 is a "good" seed for the Mersenne Twister (it is "valid", but that's not the same thing), so in RandStream, we chose to make 0 map to the one that the designers chose as their default, because people like to type 0 and we want that to be "good". Presumably you used zero because that is the "default" seed for the even older generators activated by 'seed' and 'state' (in fact, those zeros also really map to something other than zero, for the same reason). Is this a real concern, i.e., that you can't reproduce those same numbers for that one seed, or just something you noticed?
>
> Pick any other seed, or any other of the old generators, and results are identical.
>
> >> rand('twister', 5489); rand(1,3)
> ans =
> 0.81472 0.90579 0.12699
> >> h = RandStream('mt19937ar', 'Seed', 5489); rand(h, 1, 3)
> ans =
> 0.81472 0.90579 0.12699
>
> >> rand('state', 0); rand(1,3)
> ans =
> 0.95013 0.23114 0.60684
> >> h = RandStream('swb2712', 'Seed', 0); rand(h, 1, 3)
> ans =
> 0.95013 0.23114 0.60684
From: Kevin on 17 Mar 2010 01:15
Great. I think I got it now on (5).

> > You might then say, "Isn't that what seeds are for?", and I direct you to the blog entries I cited in my full response.

I guess my answer is that it is not a good idea to use different seeds to generate different random numbers (it works) but the random numbers may not have a good statistical properties. I trust that the algorithm will do the right thing and all I do is to save-and-restore the state when I generate the random numbers.

Peter Perkins <Peter.Perkins(a)MathRemoveThisWorks.com> wrote in message <hnoro9$89t$1(a)fred.mathworks.com>...
> On 3/16/2010 2:40 AM, Kevin wrote:
>
> > (5) This is about the concept of Substream in those generator that
> > supports multiple streams and substream. I admit that I know very little
> > about random number generator algorithm. I would hope that the
> > documentation would tell me the interaction between setting Substream
> > and setting State in those generators. For example, if I want to
> > generate the same random numbers in those generators (i.e. mlfg6331_64
> > and mrg32k3a), do I need to save-and-set both substream index and state
> > or just either one of them?
>
> This one I know I answered on another thread. In short, either, not both. Setting the substream index alters the state. Each substream index is in effect an alias for a specific internal state.
From: Kevin on 17 Mar 2010 11:50
Peter,

Thanks a lot for all the information being shared here. I learned a great deal. This is a very good discussion on RandStream.
From: Peter Perkins on 17 Mar 2010 13:25
On 3/17/2010 1:10 AM, Kevin wrote:

> Now I think about this more. I get a bit confused. Do both
> RandStream('mt19937ar') and rand('twister', seed) use the same Mersenne
> twister algorithm? Or it is the mapping of 0 being different in the old
> and new ways?

There are more than one Mersenne Twister algorithms (which is why we used the apparently cryptic string 'mt19937ar' to identify this ONE -- it is completely unambiguous, google it and see). But the old rand('twister') and the new RandStream('mt19937ar') are the same one, and the thing you noticed about a zero seed is the ONLY difference.

> This is a concern to me because I would like to replace the old rand and
> randn with RandStream and at the same time I want to make sure that I
> get the same random numbers.

I can't say whether or not this is an appropriate requirement, but unless you are trying to replicate a benchmark test, it seems unnecessary. Most MATLAB users _won't_ need to worry about this.

If the need to replicate exactly _is_ for a benchmark, then you have to consider that RAND and RANDN are now drawing from the same stream. In addition, while the Mersenne Twister that you get at MATLAB startup as of R2008b replicates the values from RAND exactly, the values from RANDN are different -- you'd need to use a stream created as RandStream('shr3cong') to replicate the old RANDN behavior. If you only need to replicate one or the other, or you can live with using two streams, then do that, but otherwise, consider using rand('twister',...) and randn('state',...) and ignoring the M-Lint warnings.

Again, none of the above is important unless you really care about getting exactly the same random numbers as you did prior to 2008b, and in most cases, you shouldn't.
First  |  Prev  | 
Pages: 1 2 3
Prev: How to run Matlab 6.5 in Windows XP x64 Edition
Next: nested simulation or time-generating "for"