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