We show how to generate pseudo-random $Unif[0,1]$ in R. The process is very similar in Python. The theory is discussed in Lange 22.1 and 22.2. The Python code for random number generation is here. We focus on what is going on behind the scenes in R.
Later we show how to transform $Unif[0,1]$ deviates to generate random numbers from other distributions, such as the bernoulli, exponential, and normal.
ls(all.names=TRUE)
options(repr.plot.height=4,repr.plot.width=5)
ls(all.names=TRUE)
## R makes generating random numbers look effortless
runif(1)
runif(1)
?.Random.seed
## view code of function by executing function without ()
## R calls C functions in many cases
runif
## View C code here:
## https://github.com/wch/r-source/blob/trunk/src/nmath/runif.c
function (n, min = 0, max = 1)
.Call(C_runif, n, min, max)
## .Random.seed contains RNG state
## .Random.seed[1] codes which random number generator
## first digit is generator for normal, second and third
## are for uniform, 00 = Wichmann-Hill
head(.Random.seed)
length(.Random.seed)
## usually a good idea to set the seed once (but only once)
## in a simulation
## ensures randomness and reproducibility
set.seed(1)
runif(1)
set.seed(1)
runif(1)
### default RNG is "Mersenne twister: a 623-dimensionally equidistributed
### uniform pseudo-random number generator"
### as the name suggests, it is fairly complicated
RNGkind()
### we study the Wichmann-Hill RNG instead
## much simpler but gives the main ideas of what is going on
RNGkind("Wichmann-Hill")
RNGkind()
## how does set.seed set the state in .Random.seed
.Random.seed
## when you run set.seed(some number), this sets .Random.seed to new values
set.seed(1234)
.Random.seed
runif(1)
.Random.seed
runif(1)
set.seed(1234)
.Random.seed
runif(1)
runif(1)
When running simulation studies it is a good idea to set the seed once at the beginning of the program so that your results can be reproduced exactly. If you do not set the seed manually the seed is usually chosen based on the current time.
The next "random" number is a deterministic function of the current seed values the deterministic function mimicks randomness now show how this is done for Wichmann-Hill (WH). WH uses two ideas:
p <- 30269
a <- 171
seed <- 15
gen_random_unif <- function(){
## <<- is poor coding style, don't use in practice
## this sets the value of seed *outside* of the function
seed <<- (a*seed) %% p
return(seed/p)
}
gen_random_unif()
gen_random_unif()
gen_random_unif()
## run this function 1000 times
x <- replicate(1000,gen_random_unif())
options(repr.plot.height=4,repr.plot.width=5)
hist(x)
mean(x)
High quality random number generators will pass the Diehard tests
### above mutliplicative generator repeats itself
### after p-1 random numbers
x <- replicate(p,gen_random_unif())
head(x)
tail(x)
hist(x)
## more realistic level of uncertainty
hist(runif(p))
Idea: Generate 3 random uniform random deviates using 3 different multiplicative generators (3 $a$ values and 3 $p$).
## following code taken from examples in ?.Random.seed
## This shows how 'runif(.)' works for Wichmann-Hill,
## using only R functions:
p.WH <- c(30269, 30307, 30323)
a.WH <- c( 171, 172, 170)
next.WHseed <- function(i.seed = .Random.seed[-1]) {
(a.WH * i.seed) %% p.WH
}
.Random.seed
.Random.seed[-1]
next.WHseed()
runif(n=1)
.Random.seed
## "predict" the next random number generated by R
## generates n random uniform from Wichmann-Hill
## using starting seed
## does not alter seed, so runif (built in)
## R random number generator will give same result
myunif <- function(n,seed){
p.WH <- c(30269, 30307, 30323)
a.WH <- c( 171, 172, 170)
X <- matrix(0,ncol=3,nrow=n+1)
X[1,] <- seed
for(ii in 2:(n+1)){
X[ii,] <- (a.WH * X[ii-1,]) %% p.WH
}
x <- colSums(t(X) / p.WH) %% 1
return(x[-1])
}
n <- 10
myunif(n,.Random.seed[-1])
runif(n) ## updates .Random.seed, so next run differet
### wichmann hill has a cycle length approx 6.95×10^12