quasiRNG {randtoolbox} | R Documentation |
the Torus algorithm, the Sobol and Halton sequences.
halton(n, dim = 1, init = TRUE, normal = FALSE, usetime = FALSE, mixed = FALSE, method="C", mexp = 19937) sobol(n, dim = 1, init = TRUE, scrambling = 0, seed = 4711, normal = FALSE, mixed = FALSE, method="Fortran", mexp = 19937) torus(n, dim = 1, prime, init = TRUE, mixed = FALSE, usetime = FALSE, normal=FALSE, mexp = 19937)
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 |
normal |
a logical if normal deviates are needed, default |
scrambling |
an integer value, if 1, 2 or 3 the sequence is scrambled
otherwise not. If |
seed |
an integer value, the random seed for initialization
of the scrambling process. By default |
prime |
a single prime number or a vector of prime numbers to be used in the Torus sequence. (optional argument). |
mixed |
a logical to combine the QMC algorithm with the SFMT algorithm, default |
usetime |
a logical to use the machine time to start the Torus sequence,
default |
method |
a character string either |
mexp |
an integer for the Mersenne exponent of SFMT algorithm,
only used when |
The currently available generator are given below.
Whatever the sequence, when normal=TRUE
, outputs are transformed with
the standard normal quantile (qnorm
).
The kth term of the Torus algorithm in d dimension is given by
u_k = (frac(k sqrt(p_1)), ..., frac(k sqrt(p_d)) )
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 https://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 so 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
or a randomized when mixed = TRUE
.
Computes uniform Sobol low discrepancy numbers.
The sequence starts from k=1 when initialized with
init = TRUE
(default).
When scrambling > 0
, a scrambling is performed
or when mixed = TRUE
, a randomized seed is performed.
If some number of Sobol sequences are generated outside (0,1) with scrambling,
the seed is randomized until we obtain all numbers in (0,1).
Two versions of Sobol sequences are available the current version in
Fortran (method = "Fortran"
) and the development version in C (method = "C"
).
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
(default) and no not depending on machine time
usetime = FALSE
.
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
.
Two versions of Halton sequences are available the historical version in
Fortran (method = "Fortran"
) and the new version in C (method = "C"
).
If method = "C"
, mixed
argument can be used to randomized
the Halton sequences.
See the pdf vignette for details.
torus
, halton
and sobol
generates random variables in (0,1). It returns a nxdim matrix, when dim>1
otherwise a vector of length n
.
Christophe Dutang and Diethelm Wuertz
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)
pseudoRNG
for pseudo random number generation, .Random.seed
for what is done in R about random number generation.
# (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)