quasiRNG: Toolbox for quasi random number generation
In randtoolbox: Toolbox for Pseudo and Quasi Random Number Generation and Random Generator Tests
quasiRNG R Documentation
Toolbox for quasi random number generation
Description
the Torus algorithm, the Sobol and Halton sequences.
Usage
torus(n, dim = 1, prime, init = TRUE, mixed = FALSE, usetime = FALSE,
normal = FALSE, mexp = 19937, start = 1)
sobol(n, dim = 1, init = TRUE, scrambling = 0, seed = NULL, normal = FALSE,
mixed = FALSE, method = "C", mexp = 19937, start = 1,
maxit = 10)
halton(n, dim = 1, init = TRUE, normal = FALSE, usetime = FALSE,
mixed = FALSE, method = "C", mexp = 19937, start = 1)
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 to the start value, 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 scrambling=1, Owen type type of scrambling is
applied, if scrambling=2, Faure-Tezuka type of scrambling, is
applied, and if scrambling=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 (only for sobol
with scrambling>0).
If NULL, set to 4711.
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 FALSE.
usetime
a logical to use the machine time to start the Torus sequence,
default FALSE, i.e. when usetime=FALSE the Torus sequence start from the first term.
usetime is only used when mixed=FALSE.
method
a character string either "C".
Note that mixed=TRUE is only available when method="C".
mexp
an integer for the Mersenne exponent of SFMT algorithm,
only used when mixed=TRUE.
start
an integer 0 or 1 to initiliaze the sequence, default to 1,
only used when init=TRUE.
maxit
a positive integer used to control inner loops both
for generating randomized seed and for controlling
outputs (when needed).
Details
Scrambling is temporarily disabled and will be reintroduced in a future release.
The currently available generator are given below.
Whatever the sequence, when normal=TRUE, outputs are transformed with
the quantile of the standard normal distribution qnorm.
If init=TRUE, the default, unscrambled and unmixed-SFMT quasi-random
sequences start from start.
If start != 0 and normal=FALSE,
we suggest to use 0 as recommended by Owen (2020).
One must handle the starting value (0) correctly if a quantile
function of a non-lower-bounded distribution is used.
- Torus algorithm:
-
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.
- (scrambled) Sobol sequences
-
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).
One version of Sobol sequences is available the current version in
Fortran (method = "Fortran") since method = "C" is under development.
- 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 (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.
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.
doi: 10.1145/42288.214372
Joe S., Kuo F.Y. (2003),
Remark on Algorithm 659: Implementing Sobol's Quaisrandom Seqence Generator,
ACM Transactions on Mathematical Software 29, 49–57.
doi: 10.1145/641876.641879
Joe S., Kuo F.Y. (2008),
Constructing Sobol sequences with better two-dimensional projections,
SIAM J. Sci. Comput. 30, 2635–2654,
doi: 10.1137/070709359.
Owen A.B. (2020), On dropping the first Sobol' point,
doi: 10.1007/978-3-030-98319-2_4.
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, mixed=TRUE)
## Not run:
par(mfrow = c(1,2))
acf(torus(10^6))
acf(torus(10^6, mixed=TRUE))
## End(Not run)
#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 = 500, dim = 1), main = "Uniform Halton",
xlab = "x", col = "steelblue3", border = "white")
hist(halton(n = 500, dim = 1, norm = TRUE), main = "Normal Halton",
xlab = "x", col = "steelblue3", border = "white")
# (3) Sobol sequences
#
# uniform variate
sobol(n = 10, dim = 5)
# normal variate
sobol(n = 10, dim = 5, normal = TRUE)
# some plots
hist(sobol(500, 1), main = "Uniform Sobol",
xlab = "x", col = "steelblue3", border = "white")
hist(sobol(500, 1, 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 a 2022 macbook (2017 macbook / 2007 macbook), mean of 1000 runs
#
## Not run:
# algorithm time in seconds for n=10^6
# Torus algo 0.00689 (0.012 / 0.058)
# mixed Torus algo 0.009354 (0.018 / 0.087)
# Halton sequence 0.023575 (0.180 / 0.878)
# Sobol sequence 0.010444 (0.027 / 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=TRUE)[3]) )
mean( replicate( 1000, system.time( halton(n), gcFirst=TRUE)[3]) )
mean( replicate( 1000, system.time( sobol(n), gcFirst=TRUE)[3]) )
## End(Not run)
randtoolbox documentation built on Jan. 29, 2023, 3:02 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| quasiRNG | R Documentation |
Toolbox for quasi random number generation
Description
the Torus algorithm, the Sobol and Halton sequences.
Usage
torus(n, dim = 1, prime, init = TRUE, mixed = FALSE, usetime = FALSE,
normal = FALSE, mexp = 19937, start = 1)
sobol(n, dim = 1, init = TRUE, scrambling = 0, seed = NULL, normal = FALSE,
mixed = FALSE, method = "C", mexp = 19937, start = 1,
maxit = 10)
halton(n, dim = 1, init = TRUE, normal = FALSE, usetime = FALSE,
mixed = FALSE, method = "C", mexp = 19937, start = 1)
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 |
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 (only for |
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 |
start |
an integer 0 or 1 to initiliaze the sequence, default to 1,
only used when |
maxit |
a positive integer used to control inner loops both for generating randomized seed and for controlling outputs (when needed). |
Details
Scrambling is temporarily disabled and will be reintroduced in a future release.
The currently available generator are given below.
Whatever the sequence, when normal=TRUE, outputs are transformed with
the quantile of the standard normal distribution qnorm.
If init=TRUE, the default, unscrambled and unmixed-SFMT quasi-random
sequences start from start.
If start != 0 and normal=FALSE,
we suggest to use 0 as recommended by Owen (2020).
One must handle the starting value (0) correctly if a quantile
function of a non-lower-bounded distribution is used.
- Torus algorithm:
-
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.1234will be cast toprime=7before computing Torus sequence). The Torus sequence starts from k=1 when initialized withinit = TRUEand so not depending on machine timeusetime = FALSE. This is the default. Wheninit = 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 withusetime = TRUE, which overridesinitor a randomized whenmixed = TRUE. - (scrambled) Sobol sequences
-
Computes uniform Sobol low discrepancy numbers. The sequence starts from k=1 when initialized with
init = TRUE(default). Whenscrambling > 0, a scrambling is performed or whenmixed = 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). One version of Sobol sequences is available the current version in Fortran (method = "Fortran") sincemethod = "C"is under development. - 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(default) and no not depending on machine timeusetime = FALSE. Wheninit = 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 withusetime = TRUE, which overridesinit. Two versions of Halton sequences are available the historical version in Fortran (method = "Fortran") and the new version in C (method = "C"). Ifmethod = "C",mixedargument can be used to randomized the Halton sequences.
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. doi: 10.1145/42288.214372
Joe S., Kuo F.Y. (2003), Remark on Algorithm 659: Implementing Sobol's Quaisrandom Seqence Generator, ACM Transactions on Mathematical Software 29, 49–57. doi: 10.1145/641876.641879
Joe S., Kuo F.Y. (2008), Constructing Sobol sequences with better two-dimensional projections, SIAM J. Sci. Comput. 30, 2635–2654, doi: 10.1137/070709359.
Owen A.B. (2020), On dropping the first Sobol' point, doi: 10.1007/978-3-030-98319-2_4.
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, mixed=TRUE) ## Not run: par(mfrow = c(1,2)) acf(torus(10^6)) acf(torus(10^6, mixed=TRUE)) ## End(Not run) #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 = 500, dim = 1), main = "Uniform Halton", xlab = "x", col = "steelblue3", border = "white") hist(halton(n = 500, dim = 1, norm = TRUE), main = "Normal Halton", xlab = "x", col = "steelblue3", border = "white") # (3) Sobol sequences # # uniform variate sobol(n = 10, dim = 5) # normal variate sobol(n = 10, dim = 5, normal = TRUE) # some plots hist(sobol(500, 1), main = "Uniform Sobol", xlab = "x", col = "steelblue3", border = "white") hist(sobol(500, 1, 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 a 2022 macbook (2017 macbook / 2007 macbook), mean of 1000 runs # ## Not run: # algorithm time in seconds for n=10^6 # Torus algo 0.00689 (0.012 / 0.058) # mixed Torus algo 0.009354 (0.018 / 0.087) # Halton sequence 0.023575 (0.180 / 0.878) # Sobol sequence 0.010444 (0.027 / 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=TRUE)[3]) ) mean( replicate( 1000, system.time( halton(n), gcFirst=TRUE)[3]) ) mean( replicate( 1000, system.time( sobol(n), gcFirst=TRUE)[3]) ) ## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.