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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._