RNG initialization
Paul Gilbert
pgilbert@bank-banque-canada.ca
Tue, 18 Jan 2000 10:26:26 -0500
>Since set.seed seems a suitable name and it would help compatibility,
>would
>set.seed(seed, kind = 0)
>be a suitable interface? Or would people prefer
>RNGinit(kind, seed)
Before making too many changes to random number generation it might be useful to
look at the code below which is part of my syskern package. It sets the seed,
uniform RNG and the normal transformation, and may serve as a good model for
setting other transformation. Since it also works in Splus it allows me to run
the same random experiments in Splus and R. I've been using it successfully for
some time now, on much more complicated examples than the simple tests below.
Paul Gilbert
##############################################################################
# 6/ Random number generation.
# part of syskern package. I think this is self contained except for is.S
is.S <- function(){FALSE}
##############################################################################
if (is.S())
{
RNGkind <- function(kind=NULL)
{# With a null argument this returns the current kind.
# kind = "default" resets the RNG to the Splus default.
# Splus does not allow arbitrary setting of .Random.seed so
# .RandomSeed is used. The existence of .RandomSeed is
# used to indicate whether an alternate RNG is used and the
# first element of .RandomSeed indicates the generator.
if ( is.null(kind))
{if(exists(".RandomSeed")) kind <- c("Wichmann-Hill")[1+.RandomSeed[1]]
else kind <- "default"
}
else
{# set the RNG kind
if (kind == "default")
{if(exists(".RandomSeed")) remove(".RandomSeed", where = 1) }
else if (kind == "Wichmann-Hill")
assign(".RandomSeed", c(0, as.integer(100000*runif(3))), where=1)
else stop("Only Wichmann-Hill, default or NULL supported for kind.")
}
kind
}
if (!exists("set.seed.Splus")) set.seed.Splus <- set.seed
set.seed <- function(seed=NULL)
{# with a null argument this also serves as get.seed.
kind <- RNGkind()
if ( is.null(seed))
{if (kind == "default") seed <-.Random.seed
else seed <-.RandomSeed
}
else
{# set seed
if (kind == "default")
{if (1==length(seed)) set.seed.Splus(seed)
else assign(".Random.seed", seed, where=1)
}
else if (kind == "Wichmann-Hill")
{if (3 == length(seed)) seed <- c(0,seed)
if (0 != seed[1])
stop("seed[1] should be 0 for Wichmann-Hill RNG.")
if (4 != length(seed))
stop("seed length is not consistent with kind Wichmann-Hill.")
assign(".RandomSeed", seed, where=1)
}
else stop("seed does not match RNG kind.")
}
seed
}
set.RNORMtransform <- function(rnorm=NULL)
{# with a null argument this returns setting
if (is.null(rnorm))
{if(exists(".RNORMtransform", where=1)) rnorm <- .RNORMtransform
else rnorm <- "default"
}
else
{if (rnorm == "Box-Muller") assign(".RNORMtransform", rnorm, where=1)
else if (rnorm == "default")
{if(exists(".RNORMtransform", where=1))
remove(".RNORMtransform", where=1)}
else stop("Only Box-Muller, default or NULL supported for rnorm.")
}
rnorm
}
} # end of if is.S
if (is.R())
{
# the only change in RNGkind is to allow "default" (for "Wichmann-Hill")
if (!exists("RNGkind.R")) RNGkind.R <- RNGkind
RNGkind <- function (kind = NULL)
{if ((!is.null(kind)) && (kind == "default")) kind <- "Wichmann-Hill"
RNGkind.R(kind=kind)
}
set.seed <- function(seed=NULL)
{# with a null argument this also serves as get.seed.
kind <- RNGkind()
if ( is.null(seed))
{if (!exists(".Random.seed")) z <- runif(1) # initializes seed
seed <-.Random.seed
}
else
{# set seed
if (kind == "Wichmann-Hill")
{if (3 == length(seed)) seed <- c(0,seed)
if (0 != seed[1])
stop("seed[1] should be 0 for Wichmann-Hill RNG.")
if (4 != length(seed))
stop("seed length is not consistent with kind Wichmann-Hill.")
assign(".Random.seed", seed, envir=.GlobalEnv)
}
else
{# seed <- c(kind,seed) ????
# should do some other checks or are they done elsewhere?
assign(".Random.seed", seed, envir=.GlobalEnv)
}
# else stop("seed does not match RNG kind.")
}
seed
}
set.RNORMtransform <- function(rnorm=NULL)
{# with a null argument this returns setting
if (is.null(rnorm))
{if(exists(".RNORMtransform", where=1)) rnorm <- .RNORMtransform
else rnorm <- "default"
}
else
{if (rnorm == "Box-Muller")
assign(".RNORMtransform", rnorm, envir=.GlobalEnv)
else if (rnorm == "default")
{if(exists(".RNORMtransform", where=1))
remove(".RNORMtransform", pos=1)}
else stop("Only Box-Muller, default or NULL supported for rnorm.")
}
rnorm
}
} # end of if is.R
# in both Splus and R:
set.RNG <- function(kind=NULL, seed=NULL, rnorm=NULL)
{# with a null argument this also serves as get.RNG so it is possible to
# old <- set.RNG()
# on.exit(set.RNG(old))
# set.RNG("default") can affect both kind and seed, but not rnorm
if (is.list(kind))
list(kind=RNGkind(kind$kind), seed=set.seed(kind$seed),
rnorm= set.RNORMtransform(kind$rnorm) )
else
list(kind=RNGkind(kind), seed=set.seed(seed),
rnorm=set.RNORMtransform(rnorm) )
}
if (!exists("runif.default")) runif.default <- runif
runif <- function(n, min=0, max=1, compiled=is.R())
{# This typically just calls runif.default, but allows using other
# RNGs to generate the same sequence in R and S.
# eg: set.RNG(seed=c(1:3), kind="Wichmann-Hill")
# runif(10)
if(RNGkind() == "default") return(runif.default(n, min=min, max=max))
else if(compiled) return(runif.default(n, min=min, max=max))
else seed <- set.seed() # returns the setting
kind <- c("Wichmann-Hill")[1+seed[1]]
if(kind == "Wichmann-Hill")
{out <- numeric(n)
if (seed[1] != 0) stop("seed setting is not consistent with RNG.")
x <- seed[2]; y <- seed[3]; z <- seed[4]
for(i in 1:length(out))
{x <- (171*x) %% 30269
y <- (172*y) %% 30307
z <- (170*z) %% 30323
out[i] <- (x/30269 + y/30307 + z/30323) %% 1.0
}
set.seed( c(0, x,y,z))
}
else stop("runif RNG kind not supported.")
out
}
if (!exists("rnorm.default")) rnorm.default <- rnorm
rnorm <- function(n, mean=0, sd=1, compiled=F)
{# This typically just calls rnorm.default, but provides the possibility of
# using "Wichmann-Hill" to generate the same runif sequence in R and S and
# then generate the same normally distributed numbers.
# eg: set.RNG(seed=c(1:3), kind="Wichmann-Hill", rnorm="Box-Muller")
# rnorm(10)
transf <- set.RNORMtransform()
if(transf == "default") return(rnorm.default(n, mean=mean, sd=sd))
else if(compiled) return(rnorm.default(n, mean=mean, sd=sd))
else if(transf == "Box-Muller")
{out <- runif(n + n%%2)
out <- matrix(out, 2, length(out)/2)
out <- c( sqrt(log(1/out[1,])) * sin(2*pi*out[2,]),
sqrt(log(1/out[1,])) * cos(2*pi*out[2,]))
if (1 == n%%2) out <- out[-length(out)] #drop last point if n is odd
}
else stop("normal rv transformation technique not supported.")
mean + out*sd
}
#########################################################
# test function
#########################################################
random.number.test <- function()
{cat("Random number generator tests ...")
if (is.R())
{test.seed<- c( 979, 1479, 1542) # values from 0.49 beta
test.valueU <-c(5.693354055333957e-01,1.051357751852140e-01,
5.846933178718317e-02, 7.537960906527452e-02, 7.043734921992200e-01)
test.valueN <-c(-5.559389931781886e-01,
-1.902431069568611e+00, 1.524595894866778e+00,
-7.863494805034426e-01, 1.328128164898773e-01)
}
if (is.Splus())
{test.seed<- c(37, 39, 39, 4, 7, 2, 27, 58, 38, 15, 32, 2)
test.valueU <- c(0.4299328043125570, 0.3092006836086512,
0.5808096211403608, 0.3061958812177181, 0.8137333435006440)
test.valueN <- c( -0.7613318231781665, -0.5724360196433543,
0.8536399448225964, -0.2269096022522968, -0.8126790170570223)
}
old.seed <- set.RNG(kind="default", seed=test.seed, rnorm="default")
on.exit(set.RNG(old.seed))
ok <- TRUE
if (1e-14 < max(abs(runif(5)-test.valueU)))
{warning("The default runif number generator has been changed.")
ok <- FALSE
}
set.RNG(kind="default", seed=test.seed, rnorm="default")
if (1e-14 < max(abs(rnorm(5)-test.valueN)))
{warning("The default rnorm number generator has been changed.")
ok <- FALSE
}
set.RNG(kind="Wichmann-Hill", seed=c(979,1479,1542), rnorm="Box-Muller")
if (1e-14 < max(abs(runif(5) -
c(0.56933540553339546, 0.10513577518521355, 0.05846933178718317,
0.07537960906527452, 0.70437349219921996))))
{warning("The Wichmann-Hill runif number generator has been changed.")
ok <- FALSE
}
set.RNG(kind="Wichmann-Hill", seed=c(979,1479,1542), rnorm="Box-Muller")
if (1e-14 < max(abs(rnorm(5) -
c(0.4605069059114530, 0.7685565310963474, -0.3737680932387061,
0.5926372779538560, 1.4995245125275518))))
{warning("The Box-Muller rnorm number generator has been changed.")
ok <- FALSE
}
if (ok) cat("ok\n")
else cat("failed!\n")
invisible(ok)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._