[Rd] R-devel RNG change
Paul Gilbert
pgilbert@bank-banque-canada.ca
Mon Mar 3 18:02:02 2003
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.