Re: rand("state")

[Mike Miller]

On Thu, 16 Feb 2006, Paul Kienzle wrote:

> > If you will be producing many random replicates in a simulation study, 
> > it is better to use the full 625-integer seed so that the probability 
> > of any two replicates being identical is extremely small.
> >
> > Do I have this right?
>
> All you have to worry about is the probability of your seeds colliding. 
> If you are doing 10^6 simulations, an integer seed gives you a 1 in 4000 
> chance of this happening.  If even this is too high, you can either 
> start at seed 0 and increment to get 10^6 independent sequences, or use 
> two integers for the seed, giving you a 1 in 10^13 chance of a 
> collision.
>
> You will not need a full 624 integer seed.

That helps a lot because I didn't know that consecutive integers produce 
independent sequences.  How do we know that is true?

Use of a short seed always causes the same longer state to be generated, I 
assume, and this is not dependent on machine architecture and it doesn't 
use urandom.  That will be especially helpful in saving space if many 
states have to be stored.

I don't agree with your collision probabilities.  How did you compute 
them?  This is how I would compute them:

For the probability of at least one collision in N replicates using simple 
random sampling of 2^32 seeds, we have to figure out the probability that 
there are no collisions and subtract that from 1.  Imagine doing this 
sequentially.  What is the probability that the second seed does not equal 
the first one?

(2^32-1)/2^32

If the first two are different, what is the probability that the third one 
doesn't equal the first two?

(2^32-2)/2^32

If the first N-1 integers are all different from one another, what is the 
probabiilty that the Nth one is also different from the first N-1?

(2^32-N+1)/2^32

These are independent conditional probabilities that must be multiplied 
together to get the joint probability of no collisions:

prod( ( 2^32 - [1:(N-1)] ) / 2^32 )

So the probability of at least one collision is 1 minus that probability:

1 - prod( ( 2^32 - [1:(N-1)] ) / 2^32 )

Here are some results:

octave:304> N=1000; 1-prod((2^32-[1:(N-1)])/2^32)
ans =  0.00011629215307196
octave:305> N=10000; 1-prod((2^32-[1:(N-1)])/2^32)
ans = 0.0115728899862177
octave:306> N=100000; 1-prod((2^32-[1:(N-1)])/2^32)
ans = 0.687812281963702
octave:307> N=1000000; 1-prod((2^32-[1:(N-1)])/2^32)
ans = 1

So I think the probability of at least one collision in 1,000,000 
replicates using a single 32-bit seed is approximately 1.0 (I actually get 
1-2.73901474761395e-51), and not 1/4000.  The probability of at least one 
collision exceeds 1% when there are 10,000 replicates.

This is an example of a so-called birthday problem:

http://mathworld.wolfram.com/BirthdayProblem.html

This is the answer to the original birthday problem:

octave:309> N=23; 1-prod((365-[1:(N-1)])/365)
ans = 0.507297234323986

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