R-alpha: set.seed(.) [was 'compatibility']

[Martin Maechler]

>>>>> Thomas Lumley writes:
    >> Two compatibility issues found while trying to convert a simulation
    >> from S to R.

    >> 1. set.seed()  We don't have this function.  According to Venables &
    >> Ripley it just picks a seed from a list of 1000 possibilities. How about
    >>
    >> "set.seed" <-function (i) 
    >> {
    >>         if (any(i > 1000) | any(i < 1)) 
    >>                 stop("argument to set.seed() must be in 1:1000")
    >>         if (any(i != trunc(i))) 
    >>                 stop("argument to set.seed()must be an integer")
    >>         .Random.seed <<- as.integer((200 * c(i, i^2, i^3))%%c(30269, 
    >>                 30307, 30323))
    >> }

Thank you Thomas,
yes I agree that we should have a  set.seed(.).

{{now follows a longish explanation about what S-plus does, etc..
  at the end a better (?) proposal for set.seed 
}}

When I first thought about implementing  set.seed(.) in R,
I tried to understand what S-plus does (and did not succeed).
In the mean time, I have re-investigated and found almost ``everything'' :

S-plus 3.4 *help* says

Sp3.4>>     DESCRIPTION:
Sp3.4>> 	   Puts the random number generator in a reproducible state.
Sp3.4>> 
Sp3.4>>     USAGE:
Sp3.4>> 	   set.seed(i)
Sp3.4>> 
Sp3.4>>     REQUIRED ARGUMENTS:
Sp3.4>>     i:     be an integer between 0 and 1000.
Sp3.4>> 
Sp3.4>>     SIDE EFFECTS:
Sp3.4>> 	   Sets the value of the .Random.seed object in  the  working
Sp3.4>> 	   directory.
Sp3.4>> 
Sp3.4>>     DETAILS:
Sp3.4>> 	   Random number generators in S-PLUS are all  based  upon  a
Sp3.4>> 	   single  uniform  random  number  generator  that generates
Sp3.4>> 	   numbers in a very long, but ultimately periodic  sequence.
Sp3.4>> 	   The   position   in   the   sequence  is  held  in  object
Sp3.4>> 	   .Random.seed.  Function set.seed sets .Random.seed so that
Sp3.4>> 	   subsequent  calls  to  random  number  generator functions
Sp3.4>> 	   (runif, rnorm, etc.)  will generate  numbers  from  a  new
Sp3.4>> 	   portion of the overall cycle.
Sp3.4>> 
Sp3.4>> 	   Random  number  generation  in  S-PLUS  is  adapted   from
Sp3.4>> 	   Marsaglia  (1973);  see  Kennedy  and  Gentle  (1980)  for
Sp3.4>> 	   background information.

So they __say__   ``i in 0:1000'', but the
problem is that set.seed(i) also ``works'' with any other integer i.

But what exactly should happen?
The crucial sentence is (as vague as)
 Sp3.4>> Function set.seed sets .Random.seed so that
 Sp3.4>> subsequent calls to random number generator functions (runif...) 
 Sp3.4>> will generate numbers from a new portion of the overall cycle.

The intent being that for ANY i,j (i!=j) set.seed(i) and set.seed(j) would
choose ``distant'' parts of the periodic cycle.

Unfortunately, we don't know what exactly happens in
the S-plus code for set.seed(i), namely,  .C("setseed", as.integer(i))
However, I have found out the following which should suffice :

 1) Only .Random.seed[4:6]  are varyied by set.seed, i.e., 
   after set.seed(i), .Random.seed always is

	   1  2  3  4  5  6  7  8  9 10 11 12	
      [1] 21 14 49  *  *  * 32 22 36 23 28  3
		   X4 X5 X6

   Which means that only the 'congrval' {Venables & Ripley, 2nd ed., p.166-167}
   is varied	

 2) X4 is always ==  16*((i+2) %% 4)
    X5 takes (all) values in 0:63  [also with some periodicity]
    X6 takes the values of 0:3

 3) set.seed(i) <==> set.seed(i + n*1024)  for all integers n, i

 4) { set.seed(i) ;  i \in {integers}} 
    results in a set of 1024 different values of 'congrval'.
    These values are EXACTLY EQUISPACED lying in {0,..,2^32}.

    If C[1:1024] are these sorted values, we have

	C[i] == 201621 + (i-1)* 2^22  , for i in 1:1024
	


So far for S-plus.

-----------------------------------------------------------------

Now, in R with a different basic RNG (with a longer period), 
		-- see ?.Random.seed --
we would want  set.seed(i)
to also guarantee that the locations in the periodic cycle are far apart.

I think (but may be wrong; I must do other things now !!)
this can be better achieved by just changing one of the 3 congruential
generators involved.
This leads to the following proposal :

set.seed <- function (i) 
{
        if(length(i) != 1  || i != (ii <- as.integer(i)))
	  stop("argument to set.seed() must be integer(1)")
	i <- ii
        if (i > 1024 || i < 1) {
	  warning("set.seed(i) has i outside 1:1024; taking  i %% 1024")
	  i <- i %% 1024; if(i <= 0) i <- i+ 1024
	}
        invisible(.Random.seed <<- as.integer(c((171*i)%%30269, 10909, 10111)))
}



Martin Maechler <maechler@stat.math.ethz.ch>		 <><
Seminar fuer Statistik, SOL G1
ETH (Federal Inst. Technology)	8092 Zurich	 SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1086
http://www.stat.math.ethz.ch/~maechler/
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-