[R] random sequence
Prof Brian D Ripley
ripley at stats.ox.ac.uk
Tue Apr 27 08:53:33 CEST 1999
On Mon, 26 Apr 1999, Paul Gilbert wrote:
> Brian
>
> >What exactly do you want? My class exercise is to implement the
> >Wichmann-Hill algorithm in C/Fortran and use it with S-PLUS! But
> >while one is at it, why not use a decent generator with both?
>
> It would certainly be nice to have all the R generators in S, but for now I
> would be happy to have just one. What I would like is to be able to run the
> programs which test my code, and get the same random sequences so that I don't
> need different check values for R and Splus.
>From my files (from before R was thought of):
WichHill <- function(n)
{
output <- numeric(n)
seed <- get(".WHseed", frame=0)
x <- seed[1]; y <- seed[2]; z <- seed[3]
for(i in 1:n){
x <- (171*x) %% 30269
y <- (172*y) %% 30307
z <- (170*z) %% 30323
output[i] <- (x/30269 + y/30307 + z/30323) %% 1.0
}
assign(".WHseed", c(x,y,z), frame=0, immediate=T)
output
}
> assign(".WHseed", 1:3, frame=0)
> WichHill(10)
[1] 0.03381877 0.77754189 0.05273525 0.74462407 0.49036219 0.98285437
[7] 0.80915099 0.71338138 0.80102091 0.98958603
> WichHill(10)
[1] 0.9168563 0.4666467 0.6701995 0.9274966 0.8431496 0.7882993 0.7696486
[8] 0.2939858 0.6187752 0.9892436
in R:
> .Random.seed <- c(0, 1:3) # select Wichmann-Hill, seed = 1:3
> runif(10)
[1] 0.03381877 0.77754189 0.05273525 0.74462407 0.49036219 0.98285437
[7] 0.80915099 0.71338138 0.80102091 0.98958603
> runif(10)
[1] 0.9168563 0.4666467 0.6701995 0.9274966 0.8431496 0.7882993 0.7696486
[8] 0.2939858 0.6187752 0.9892436
You can do this faster in C/Fortran, of course, but the load/save of the
seed is slow from compiled code and in pure S I get over 1000 unifs/sec.
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help 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-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list