pseudoRNG: 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.

The list of supported generators consists of generators available via direct functions and generators available via set.generator() and runif() interface. Most of the generators belong to both these groups, but some generators are only available directly (SFMT) and some only via runif() interface (Mersenne Twister 2002). This help page describes the list of all the supported generators and the functions for the direct access to those, which are available in this way. See set.generator() for the generators available via runif() interface.

Usage

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, version = "a")
knuthTAOCP(n, dim = 1)
setSeed(seed)

Arguments

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

dim

dimension of observations (must be

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.

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.

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

version

a character either 'a' or 'b'. see details

Value

SFMT, WELL, congruRand and knuthTAOCP generate random variables in ]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 sets the seed of the randtoolbox package (i.e. both for the knuthTAOCP, 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],[object Object]See the pdf 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)

Wikipedia (2008), a linear congruential generator.

Examples

Run this code

require(rngWELL)

# (1) 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)

# (2) 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)


# (3) 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] )

# (4) 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
setSeed4WELL(08082008)
WELL(10, order =19937)
setSeed4WELL(08082008)
WELL(10, order =19937, temper=TRUE)

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


# (6) How to set the seed?
# all example is duplicated to ensure setSeed works

# congruRand
setSeed(1302)
congruRand(1)
setSeed(1302)
congruRand(1)
# SFMT
setSeed(1302)
SFMT(1, usepset=FALSE)
setSeed(1302)
SFMT(1, usepset=FALSE)
# BEWARE if you do not set usepset to FALSE
setSeed(1302)
SFMT(1)
setSeed(1302)
SFMT(1)
# WELL
setSeed(1302)
WELL(1)
setSeed(1302)
WELL(1)
# Knuth TAOCP
setSeed(1302)
knuthTAOCP(1)
setSeed(1302)
knuthTAOCP(1)



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

# algorithm			time in seconds for n=10^6
# classical Mersenne Twister  			0.066 
# SF Mersenne Twister  	       			0.044 
# WELL generator				0.065
# Knuth TAOCP					0.046
# Park Miller  		      			0.108
n <- 1e+06
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]) )

Run the code above in your browser using DataLab