From pgilbert@bank-banque-canada.ca  Mon Mar  3 18:02:02 2003
From: pgilbert@bank-banque-canada.ca (Paul Gilbert)
Date: Mon Mar  3 18:02:02 2003
Subject: [Rd] R-devel RNG change
Message-ID: <3E6371D3.4F13F1CB@bank-banque-canada.ca>

I find the documention for RNGversion in R-devel is a bit misleading, 
and suggest adding a sentence to make it clear that the meaning of 
"default" is not set to its meaning in the earlier R version:

      `RNGversion' can be used to set the random generators as they were
      in an earlier R{} version (for reproducibility). RNGversion does
      not set the meaning of "default" to its meaning in the earlier R{}
      version.

I also suggest the addition of a small function to query and set all 
aspects of random number generation. I find it is very useful to program with a
construct like:

 simulate <- function(model, rng=NULL, ...) 
   {if(is.null(rng)) rng <- set.RNG() # current settings to return with
simualtion object
    else        {old.rng <- set.RNG(rng);  on.exit(set.RNG(old.rng))  }
    ...
   }

and know that I have all the information for reproducibility. Below is a version
of this that I have been using for a few years. Its main shortcoming is that it
only considers the uniform generator and normal transformation.

Paul Gilbert
______
set.RNG <- function(kind=NULL, seed=NULL, normal.kind=NULL)
       {# with a null argument this also serves as get.RNG
        #The next line means that set.RNG with null args does not truly
        #  return the state of the RNG in the case when it has not been
        #  initialize. It first initializes it and then returns the state. The
        #  rational is that querying the state is usually for the purpose of
        #  reproducing it, so it must first be initialized to put it in a
        #  reproducible state.
        if (!exists(".Random.seed")) z <- runif(1)
        old <- list(kind=RNGkind()[1],normal.kind=RNGkind()[2], 
                    seed=.Random.seed[-1])
        if (is.null(kind) & is.null(seed) & is.null(normal.kind)) return (old)
        #kind can be a list as returned by a previous call to set.RNG
        if (is.list(kind))
           {seed        <- kind$seed
            normal.kind <- kind$normal.kind
            kind        <- kind$kind
           }
        remove(".Random.seed", envir=.GlobalEnv) # otherwise RNGkind complains
        RNGkind(kind=kind, normal.kind=normal.kind)
        if ( 1==length(seed)) set.seed(seed)
         else assign(".Random.seed", c(.Random.seed[1], seed), 
envir=.GlobalEnv)
        old
       }

# extract the RNG information from a (simulation) object
get.RNG <- function(e=NULL)UseMethod("get.RNG")
get.RNG.default <- function(e=NULL)
   {if (is.null(e)) return(set.RNG())
    # not sure if version information is useful?
    if      (!is.null(e$version))  v <- e$version
    else if (!is.null(e$noise))    v <- e$noise$version
    else v <- NULL
    if (is.null(v)) warning("version cannot be verified getting random 
seed.")
    else if (!all(unlist(version) == unlist(v)))
      warning("Seed used but version does not correspond to original. 
Differences may occur.")
   if      (!is.null(e$rng))  k <- e$rng
   else if (!is.null(e$noise)) k <- e$noise$rng
   else stop("RNG info not found.")
   k
  }

I also have documentation, but it needs to be cleaned up a bit.



[Rd] R-devel RNG change
Paul Gilbert