the Torus algorithm, the Sobol and Halton sequences.
Usage
halton(n, dim = 1, init = TRUE, normal = FALSE, usetime = FALSE)
sobol(n, dim = 1, init = TRUE, scrambling = 0, seed = 4711, normal = FALSE)
torus(n, dim = 1, prime, init = TRUE, mixed = FALSE, usetime = FALSE,
normal=FALSE)
Arguments
n
number of observations. If length(n) > 1,
the length is taken to be the required number.
dim
dimension of observations default 1.
init
a logical, if TRUE the sequence is initialized and
restarts, otherwise not. By default TRUE.
normal
a logical if normal deviates are needed, default FALSE
scrambling
an integer value, if 1, 2 or 3 the sequence is scrambled
otherwise not. If 1, Owen type type of scrambling is
applied, if 2, Faure-Tezuka type of scrambling, is
applied, and if 3, both Owen+Faure-Tezuka type of
scrambling is applied. By default 0.
seed
an integer value, the random seed for initialization
of the scrambling process. By default 4711. On effective
if scrambling>0.
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.
Details
The currently available generator are given below.
Torus algorithm:
The kth term of the Torus algorithm in d dimension is given by
where p_i denotes the ith prime number, frac the fractional part
(i.e. frac(x) = x-floor(x)). We use the 100 000 first prime numbers
from http://primes.utm.edu/, thus the dimension is limited to 100 000.
If the user supplys prime numbers through
the argument prime, we do NOT check for primality and we cast numerics
to integers, (i.e. prime=7.1234 will be cast to prime=7 before
computing Torus sequence).
The Torus sequence starts from k=1 when initialized with
init = TRUE and no not depending on machine time
usetime = FALSE. This is the default. When init = FALSE, the sequence
is not initialized (to 1) and starts from the last term. We can also use the
machine time to start the sequence with usetime = TRUE, which overrides
init.
scrambled Sobol sequences
Calculates a matrix of uniform and normal deviated Sobol low
discrepancy numbers. Optional scrambling of the sequence can
be selected.
The Sobol sequence restarts and is initialized when
init = TRUE and otherwise not.
Halton sequences
Calculates a matrix of uniform or normal deviated halton low
discrepancy numbers. Let us note that Halton sequence in dimension is the
Van Der Corput sequence.
The Halton sequence starts from k=1 when initialized with
init = TRUE and no not depending on machine time
usetime = FALSE. This is the default. When init = FALSE, the sequence
is not initialized (to 1) and starts from the last term. We can also use the
machine time to start the sequence with usetime = TRUE, which overrides
init.
See the pdf vignette for details.
Value
torus, halton and sobol generates random variables in ]0,1[. It returns a nxdim matrix, when dim>1 otherwise a vector of length n.
Author(s)
Christophe Dutang and Diethelm Wuertz
References
Bratley P., Fox B.L. (1988);
Algorithm 659: Implementing Sobol's Quasirandom
Sequence Generator,
ACM Transactions on Mathematical Software 14, 88–100.
Joe S., Kuo F.Y. (1998);
Remark on Algorithm 659: Implementing Sobol's Quaisrandom
Seqence Generator.
Planchet F., Jacquemin J. (2003), L'utilisation de methodes de
simulation en assurance. Bulletin Francais d'Actuariat, vol. 6, 11, 3-69. (available online)
See Also
pseudoRNG for pseudo random number generation, .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))
#usage of the init argument
torus(5)
torus(5, init=FALSE)
#should be equal to the combination of the two
#previous call
torus(10)
# (2) Halton sequences
#
# uniform variate
halton(n = 10, dim = 5)
# normal variate
halton(n = 10, dim = 5, normal = TRUE)
#usage of the init argument
halton(5)
halton(5, init=FALSE)
#should be equal to the combination of the two
#previous call
halton(10)
# some plots
par(mfrow = c(2, 2), cex = 0.75)
hist(halton(n = 5000, dim = 1), main = "Uniform Halton",
xlab = "x", col = "steelblue3", border = "white")
hist(halton(n = 5000, dim = 1, norm = TRUE), main = "Normal Halton",
xlab = "x", col = "steelblue3", border = "white")
# (3) Sobol sequences
#
# uniform variate
sobol(n = 10, dim = 5, scrambling = 3)
# normal variate
sobol(n = 10, dim = 5, scrambling = 3, normal = TRUE)
# some plots
hist(sobol(5000, 1, scrambling = 2), main = "Uniform Sobol",
xlab = "x", col = "steelblue3", border = "white")
hist(sobol(5000, 1, scrambling = 2, normal = TRUE), main = "Normal Sobol",
xlab = "x", col = "steelblue3", border = "white")
#usage of the init argument
sobol(5)
sobol(5, init=FALSE)
#should be equal to the combination of the two
#previous call
sobol(10)
# (4) computation times on my macbook, mean of 1000 runs
#
## Not run:
# algorithm time in seconds for n=10^6
# Torus algo 0.058
# mixed Torus algo 0.087
# Halton sequence 0.878
# Sobol sequence 0.214
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( halton(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( sobol(n), gcFirst=TRUE)[3]) )
## End(Not run)