qrng {qrng} | R Documentation |
Computing Korobov, generalize Halton and Sobol' quasi-random sequences.
korobov(n, d = 1, generator, randomize = c("none", "shift")) ghalton(n, d = 1, method = c("generalized", "halton")) sobol (n, d = 1, randomize = c("none", "digital.shift", "Owen", "Faure.Tezuka", "Owen.Faure.Tezuka"), seed, skip = 0, ...)
n |
Number n of points to be generated >= 2. |
d |
Dimension d. |
generator |
A |
randomize |
A
|
method |
A |
seed |
if provided, an integer used within
|
skip |
number of initial points in the sequence to be skipped
( |
... |
additional arguments passed to |
Note that these procedures call fast C code. The following restrictions apply:
n
,d
must be <= 2^31-1.
n
must be <= 2^32-1 and
d
must be <= 360.
if randomize = "none"
or randomize =
"digital.shift"
, n
must be <= 2^31-1 and
d
must be <= 16510.
The choice of parameters for korobov()
is crucial for the quality of
this quasi-random sequence (only basic sanity checks are
conducted). For more details, see l'Ecuyer and Lemieux (2000).
The generalized Halton sequence uses the scrambling factors of Faure and Lemieux (2009).
See the example below on being careful about using skip > 0
when randomize = TRUE
; in this case, choosing a wrong seed
(or no seed) might lead to a bad sequence.
korobov()
and ghalton()
return an
(n,d)-matrix
; for d=1 an n-vector
is returned.
Marius Hofert and Christiane Lemieux
Faure, H., Lemieux, C. (2009). Generalized Halton Sequences in 2008: A Comparative Study. ACM-TOMACS 19(4), Article 15.
l'Ecuyer, P., Lemieux, C. (2000). Variance Reduction via Lattice Rules. Stochastic Models and Simulation, 1214–1235.
Lemieux, C., Cieslak, M., Luttmer, K. (2004). RandQMC User's guide. See https://www.math.uwaterloo.ca/~clemieux/randqmc/guide.pdf
n <- 1021 # prime d <- 4 # dimension ## Korobov's sequence generator <- 76 # see l'Ecuyer and Lemieux u <- korobov(n, d = d, generator = generator) pairs(u, gap = 0, pch = ".", labels = as.expression( sapply(1:d, function(j) bquote(italic(u[.(j)]))))) ## Randomized Korobov's sequence set.seed(271) u <- korobov(n, d = d, generator = generator, randomize = "shift") pairs(u, gap = 0, pch = ".", labels = as.expression( sapply(1:d, function(j) bquote(italic(u[.(j)]))))) ## Generalized Halton sequence (randomized by definition) set.seed(271) u <- ghalton(n, d) pairs(u, gap = 0, pch = ".", labels = as.expression( sapply(1:d, function(j) bquote(italic(u[.(j)]))))) ## Randomized Sobol' sequence (with digital shift) set.seed(271) u <- sobol(n, d, randomize = "digital.shift") pairs(u, gap = 0, pch = ".", labels = as.expression( sapply(1:d, function(j) bquote(italic(u[.(j)]))))) ## Randomized Sobol' sequence (with Owen scrambling) ## Note: - This is calling randtoolbox::sobol(, scrambling = 1, seed = ...), ## - The seed needs to be passed on directly in order to be respected ## (a global set.seed() isn't) u <- sobol(n, d, randomize = "Owen", seed = 271) pairs(u, gap = 0, pch = ".", labels = as.expression( sapply(1:d, function(j) bquote(italic(u[.(j)]))))) ## Check whether a Sobol' sequence of size 2*n equals one of size n ## and, concatenated, one of size n with the first n numbers skipped f <- function(n) { set.seed(271) a <- sobol(2*n, randomize = "digital.shift") set.seed(271) b1 <- sobol(n, randomize = "digital.shift") set.seed(271) b2 <- sobol(n, randomize = "digital.shift", skip = n) all(all.equal(apply(cbind(a, c(b1, b2)), 1, diff), rep(0, 2*n))) } stopifnot(sapply(1:10, f)) # check for n = 1, ..., 10 ## Careful when using skip > 0 and randomize = TRUE => seed matters! ## Drawing all points at once (works, of course) n <- 32 set.seed(5) U.2n <- sobol(2*n, d = 2, randomize = "digital.shift") plot(U.2n, main = "All points generated at once", xlab = expression(U[1]), ylab = expression(U[2])) ## Drawing successively with the same seed (typically the right approach) set.seed(5) U.n.1 <- sobol(n, d = 2, randomize = "digital.shift") set.seed(5) # => same seed U.n.2 <- sobol(n, d = 2, randomize = "digital.shift", skip = n) U.2n.same.seed <- rbind(U.n.1, U.n.2) plot(U.2n.same.seed, main = "Drawing successively, with the same seed", xlab = expression(U[1]), ylab = expression(U[2])) stopifnot(all.equal(U.2n, U.2n.same.seed)) # sanity check ## Drawing successively but with different seeds (typically the wrong approach) set.seed(5) U.n.1 <- sobol(n, d = 2, randomize = "digital.shift", skip = 0) set.seed(22) U.n.2 <- sobol(n, d = 2, randomize = "digital.shift", skip = n) U.2n.different.seed <- rbind(U.n.1, U.n.2) plot(U.2n.different.seed, main = "Drawing successively, with different seeds", xlab = expression(U[1]), ylab = expression(U[2]))