quasiRNG | R Documentation |
the Torus algorithm, the Sobol and Halton sequences.
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)
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). |
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.
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).
One version of Sobol sequences is available the current version in
Fortran (method = "Fortran"
) since method = "C"
is under development.
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. 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)
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, 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.