randtoolbox (version 1.04)

randtoolbox: Toolbox for pseudo and quasi random number generation

Description

General linear congruential generators such as Park Miller sequence, generalized feedback shift register such as SF-Mersenne Twister algorithm and WELL generator; and a quasi random generator (pseudo random generators) and the Torus algorithm (quasi random generation).

Usage

torus(n, dim = 1, prime, mixed = FALSE, usetime = FALSE)
congruRand(n, dim = 1, mod = 2^31-1, mult = 16807, incr = 0, echo)
SFMT(n, dim = 1, mexp = 19937, usepset = TRUE, withtorus = FALSE, usetime = FALSE)
WELL(n, dim = 1, order = 512, temper = FALSE)
knuthTAOCP(n, dim = 1)
setSeed(seed)

Arguments

n

number of observations. If length(n) > 1, the length is taken to be the number required.

dim

dimension of observations (must be

prime

a single prime number or a vector of prime numbers to be used in the Torus sequence. (optional argument).

mixed

a logical to use the mixed Torus algorithm, default FALSE.

usetime

a logical to use the machine time to start the Torus sequence, default TRUE. if FALSE, the Torus sequence start from the first term.

seed

a single value, interpreted as a positive integer for the seed. e.g. append your day, your month and your year of birth.

mod

an integer defining the modulus of the linear congruential generator.

mult

an integer defining the multiplier of the linear congruential generator.

incr

an integer defining the increment of the linear congruential generator.

echo

a logical to plot the seed while computing the sequence.

mexp

an integer for the mersenne exponent of SFMT algorithm. see details

withtorus

a numeric in ]0,1] defining the proportion of the torus sequence appended to the SFMT sequence; or a logical equals to FALSE (default).

usepset

a logical to use a set of 12 parameters set for SFMT. default TRUE.

order

a positive integer for the order of the characteristic polynomial. see details

temper

a logical if you want to do a tempering stage. see details

Value

• torus, SFMT, WELL and congruRand generates random variables in ]0,1[, ]0,1[, [0,1[ and [0,1[ respectively. It returns a $n$x$dim$ matrix, when dim>1 otherwise a vector of length n.

  setSeed set the seed of the randtoolbox package (i.e. both for the torus, SFMT, WELL and congruRand functions).

Details

The currently available generator are given below. [object Object],[object Object],[object Object],[object Object],[object Object],[object Object]See the vignette for details.

References

Knuth D. (1997), The Art of Computer Programming V2 Seminumerical Algorithms, Third Edition, Massachusetts: Addison-Wesley.

Matsumoto M. and Nishimura T. (1998), Dynamic Creation of Pseudorandom Number Generators, Monte Carlo and Quasi-Monte Carlo Methods, Springer, pp 56--69. (available online)

Matsumoto M., Saito M. (2008), SIMD-oriented Fast Mersenne Twister: a 128-bit Pseudorandom Number Generator. (available online)

Paneton F., L'Ecuyer P. and Matsumoto M. (2006), Improved Long-Period Generators Based on Linear Recurrences Modulo 2, ACM Transactions on Mathematical Software. (preprint available online)

Park S. K., Miller K. W. (1988), Random number generators: good ones are hard to find. Association for Computing Machinery, vol. 31, 10, pp 1192-2001. (available online)

Planchet F., Jacquemin J. (2003), L'utilisation de methodes de simulation en assurance. Bulletin Francais d'Actuariat, vol. 6, 11, 3-69. (available online)

Wikipedia (2008), a linear congruential generator.

See Also

.Random.seed for what is done in R about random number generation.

Examples

# (1) the Torus algorithm
#
torus(100)

# example of setting the seed
setSeed(1)
torus(5)
setSeed(6)
torus(5)
#the same
setSeed(1)
torus(10)

#no use of the machine time
torus(10, use=FALSE)

#Kolmogorov Smirnov test
#KS statistic should be around 0.0019
ks.test(torus(1000), punif) 
	
#KS statistic should be around 0.0003
ks.test(torus(10000), punif) 

#the mixed Torus sequence
torus(10, mix=TRUE)
par(mfrow = c(1,2))
acf(torus(10^6))
acf(torus(10^6, mix=TRUE))

# (2) the Park Miller sequence
#
# Park Miller sequence, i.e. mod = 2^31-1, mult = 16807, incr=0
# the first 10 seeds used in Park Miller sequence
# 16807          1
# 282475249          2
# 1622650073          3
# 984943658          4
# 1144108930          5
# 470211272          6
# 101027544          7
# 1457850878          8
# 1458777923          9
# 2007237709         10
setSeed(1)
congruRand(10, echo=TRUE)

# the 9998+ th terms 
# 925166085       9998
# 1484786315       9999
# 1043618065      10000
# 1589873406      10001
# 2010798668      10002
setSeed(1614852353) #seed for the 9997th term
congruRand(5, echo=TRUE)

# (3) the SF Mersenne Twister algorithm
SFMT(1000)

#Kolmogorov Smirnov test
#KS statistic should be around 0.037
ks.test(SFMT(1000), punif) 
	
#KS statistic should be around 0.0076
ks.test(SFMT(10000), punif) 

#different mersenne exponent with a fixed parameter set
#
SFMT(10, mexp = 607, usepset = FALSE)
SFMT(10, mexp = 1279, usepset = FALSE)
SFMT(10, mexp = 2281, usepset = FALSE)
SFMT(10, mexp = 4253, usepset = FALSE)
SFMT(10, mexp = 11213, usepset = FALSE)
SFMT(10, mexp = 19937, usepset = FALSE)
SFMT(10, mexp = 44497, usepset = FALSE)
SFMT(10, mexp = 86243, usepset = FALSE)
SFMT(10, mexp = 132049, usepset = FALSE)
SFMT(10, mexp = 216091, usepset = FALSE)

#use different sets of parameters [default when possible]
#
for(i in 1:7) print(SFMT(1, mexp = 607))
for(i in 1:7) print(SFMT(1, mexp = 2281))
for(i in 1:7) print(SFMT(1, mexp = 4253))
for(i in 1:7) print(SFMT(1, mexp = 11213))
for(i in 1:7) print(SFMT(1, mexp = 19937))

#use a fixed set and a fixed seed
#should be the same output
setSeed(08082008)
SFMT(1, usepset = FALSE)
setSeed(08082008)
SFMT(1, usepset = FALSE)


# (4) withtorus argument
# 

# one third of outputs comes from Torus algorithm
u <- SFMT(1000, with=1/3)
# the third term of the following code is the first term of torus sequence
print(u[666:670] )

# (5) WELL generator
#

# 'basic' calls
# WELL512
WELL(10, order = 512)
# WELL1024
WELL(10, order = 1024)
# WELL19937
WELL(10, order = 19937)
# WELL44497
WELL(10, order = 44497)
# WELL19937 with tempering 
WELL(10, order = 19937, temper = TRUE)
# WELL44497 with tempering
WELL(10, order = 44497, temper = TRUE)

# tempering vs no tempering
setSeed(08082008)
WELL(10, order =19937)
setSeed(08082008)
WELL(10, order =19937, temper=TRUE)

# (6) Knuth TAOCP generator
#
knuthTAOCP(10)
knuthTAOCP(10, 2) 

# (7) computation times on my macbook, mean of 1000 runs
#

# algorithm			time in seconds for n=10^6
# Torus algo					0.058 
# mixed Torus algo 	       			0.087 
# classical Mersenne Twister  			0.066 
# SF Mersenne Twister  	       			0.044 
# WELL generator				0.065
# Knuth TAOCP					0.046
# Park Miller  		      			0.174 
n <- 1e+06
mean( replicate( 1000, system.time( torus(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( torus(n, mixed=TRUE), gcFirst=T)[3]) )
mean( replicate( 1000, system.time( runif(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( SFMT(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( WELL(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( knuthTAOCP(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( congruRand(n), gcFirst=TRUE)[3]) )