MATLAB Central Newsreader - recovering the jth state of the random number generator

recovering the jth state of the random number generator

Tudor — Fri, 05 Sep 2008 23:57:01 -0400

Hi,

i am trying to produce a stimulus sequence using the following procedure:

rand ('state', 0)

for i = 1:90000
randomnoisematrix = rand (45,45)
end

Once this stimulus is presented, i have to reverse correlate the response to various epochs of the stimulus. Say, to the states of the stimulus between the 43000 and 43060 iterations.

I find that if i run the iteration above repeatedly, to any given point ( say 43001) i get back the same randomnoise matrix as expected, however if i ask for
rand ('state', 43001)
randmonoisematrix = rand (45, 45) i get a completely different output.
I other words, it appears that the state 0 is not followed by the state 1, 2, 3 etc.
Does anybody know what is going on ?

The reason i want to do this is because i do not want to have to run each time through the whole sequence to get to my epoch of interest, but rather would like to "skip" ahead to the state i want (inbetween 43000 and 43060 for example).

Thanks a lot.

Tudor

Re: recovering the jth state of the random number generator

Walter Roberson — Sat, 06 Sep 2008 01:00:28 -0400

Tudor wrote:

> I find that if i run the iteration above repeatedly, to any given point
> ( say 43001) i get back the same randomnoise matrix as expected, however if i ask for
> rand ('state', 43001)
> randmonoisematrix = rand (45, 45) i get a completely different output.
> I other words, it appears that the state 0 is not followed by the state 1, 2, 3 etc.

Correct. The state is used to initialize the internal binary data used in
the pseudo-random number generator. The state is not a counter of how many
values have been produced.

> The reason i want to do this is because i do not want to have to run each time
> through the whole sequence to get to my epoch of interest, but rather would like
> to "skip" ahead to the state i want (inbetween 43000 and 43060 for example).

You will have to create your own random number generator then. For example, you
could MD5 the state number, extract some bits from the MD5 and typecast() the
result to a floating point number.

Or you could research one of the linear congruential pseudo random number generators.
Although they are usually processed sequentially, if you go back to the theory
of them, you could go directly to any particular state without a lot more work
than raising a large number to a power (the state number) modulo a large number.
You would have to use indefinite precision arithmetic to do it the naive way as
the values get *huge*, but fortunately you can take advantage of the fact that
the number and the modulus will be relatively prime (for any PRNG of maximum
length) and there are theorems that in that case allow the work to be done to
be reduced enormously (to a practical level, if you have the right software.)

--
Q = quotation(rand);
if isempty(Q); error('Quotation server filesystem problems')
else sprintf('%s',Q), end

Re: recovering the jth state of the random number generator

Tudor — Sat, 06 Sep 2008 19:59:02 -0400

Thanks a lot for your explanation. That was a crucial thing to understand. I am a biologist trying to design some physiology experiments, so my math is far from developed.

How about if i 'seed' the rand function at each iteration ?:

for i = 1:90000

rand ('state', i);
randmonoisematrix = rand (45, 45);

end

Intuitively, this the same as sampling a random space of values, and perchance the subsample would be random as well). The goal for me is not necessarily to alter the stimulus on each presentation, but rather to ensure that the "surface" of the matrix is as "diverse" as possible, in the x,y dimension, as well as in the temporal one, for the duration of the stimulus. I imagine this would not work if there is some periodicity of the sequence that can add up to the integer value.
Thanks again for your insight. Is there a refference that describes the math behind the rand and randn functions ?

Tudor

>
> > I find that if i run the iteration above repeatedly, to any given point
> > ( say 43001) i get back the same randomnoise matrix as expected, however if i ask for
> > rand ('state', 43001)
> > randmonoisematrix = rand (45, 45) i get a completely different output.
> > I other words, it appears that the state 0 is not followed by the state 1, 2, 3 etc.
>
> Correct. The state is used to initialize the internal binary data used in
> the pseudo-random number generator. The state is not a counter of how many
> values have been produced.
>
> > The reason i want to do this is because i do not want to have to run each time
> > through the whole sequence to get to my epoch of interest, but rather would like
> > to "skip" ahead to the state i want (inbetween 43000 and 43060 for example).
>
> You will have to create your own random number generator then. For example, you
> could MD5 the state number, extract some bits from the MD5 and typecast() the
> result to a floating point number.
>
> Or you could research one of the linear congruential pseudo random number generators.
> Although they are usually processed sequentially, if you go back to the theory
> of them, you could go directly to any particular state without a lot more work
> than raising a large number to a power (the state number) modulo a large number.
> You would have to use indefinite precision arithmetic to do it the naive way as
> the values get *huge*, but fortunately you can take advantage of the fact that
> the number and the modulus will be relatively prime (for any PRNG of maximum
> length) and there are theorems that in that case allow the work to be done to
> be reduced enormously (to a practical level, if you have the right software.)

Re: recovering the jth state of the random number generator

Peter Perkins — Tue, 09 Sep 2008 14:30:40 -0400

Tudor wrote:

> How about if i 'seed' the rand function at each iteration ?:
>
> for i = 1:90000
>
> rand ('state', i);
> randmonoisematrix = rand (45, 45);
>
> end

Tudor, I recommend not doing that. Some people try to do that each time they generate a single value, and that's definitely a bad idea, because the statistical properties of the generator are completely destroyed. Reseeding each time you generate a 45x45 matrix isn't quite as bad, but still ...

What you're looking for is often called substreams.

If you have access to the MATLAB prerelease, I recommend looking at that.

Re: recovering the jth state of the random number generator

Tudor — Tue, 09 Sep 2008 16:14:02 -0400

Hi Peter,

i studied a little into it, and begin to understand how the random number generator actually works. Since it is not simply a f(x), with x being the successive states, one cannot recover the jth state after a given beginning seed except by running the actual generator with the correct seed. Since my application is to produce a succesion of matrices (which will become images after a lookup table application) that are random with respect to each other, but easy to recover to be able to run reverse correlations of stimulus responses on the given stimulus
  rand('state', myseed);
  for i = 1:framenumber
  noisematrix = rand (side,side)
  end

i have take the typical approach of the molecular biologist.
I ran the above loop offline, and stored all the frames. Since i need a grayscale lookup table (0-255), that fits in a uint8, the whole thing is under 200 Mb. Then i can reverse correlate at ease, on the saved file.

I have to say everybody that answered to my quest was very helpful, and i will be back for sure with more simpleminded questions.

Tudor

Peter Perkins <Peter.PerkinsRemoveThis@mathworks.com> wrote in message <ga61ag$1lt$1@fred.mathworks.com>...
> Tudor wrote:
>
> > How about if i 'seed' the rand function at each iteration ?:
> >
> > for i = 1:90000
> >
> > rand ('state', i);
> > randmonoisematrix = rand (45, 45);
> >
> > end
>
> Tudor, I recommend not doing that. Some people try to do that each time they generate a single value, and that's definitely a bad idea, because the statistical properties of the generator are completely destroyed. Reseeding each time you generate a 45x45 matrix isn't quite as bad, but still ...
>
> What you're looking for is often called substreams.
>
> If you have access to the MATLAB prerelease, I recommend looking at that.