RNG R/Splus compatibility
Paul Gilbert
pgilbert@bank-banque-canada.ca
Wed, 05 May 1999 12:31:59 -0400
Starting with example Wichmann-Hill code from Brian Ripley I have been playing
with a set of programs for getting the same random sequences from R and Splus. A
copy is included below along with a test (which works in Solaris with R and
Splus 3.3).
The approach is somewhat different from the usual problems on this list as I am
trying to get the same results from Splus as I get from R. However, there are a
couple of R issues:
- There is a small change to RNGkind to allow "default" (currently
Wichmann-Hill but could change from time to time).
-The strategy for specifying the transformation to get rnorm results should be
examined, as I think it should be put into R and extended for other random
number transformations.
Paul Gilbert
________
# First define these 3 functions appropriately if you don't have them already:
# is.S <- function(){T or F}
# is.R <- function(){T or F}
# is.Splus <- function(){T or F}
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))
{kind <- RNGkind(kind$kind)
seed <- set.seed(kind$seed)
rnorm <- set.RNORMtransform(kind$rnorm)
}
else
{kind <- RNGkind(kind)
seed <- set.seed(seed)
rnorm <- set.RNORMtransform(rnorm)
}
list(kind=kind, seed=seed, rnorm=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.seed(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.seed(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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._