Title: | (Randomized) Quasi-Random Number Generators |
---|---|
Description: | Functionality for generating (randomized) quasi-random numbers in high dimensions. |
Authors: | Marius Hofert [aut, cre], Christiane Lemieux [aut] |
Maintainer: | Marius Hofert <[email protected]> |
License: | GPL-2 | GPL-3 |
Version: | 0.0-10 |
Built: | 2024-11-23 05:18:18 UTC |
Source: | https://github.com/cran/qrng |
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, ...)
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 |
d |
dimension |
generator |
|
randomize |
|
method |
|
seed |
if provided, an integer used within
|
skip |
number of initial points in the sequence to be skipped
( |
... |
additional arguments passed to |
For sobol()
examples see demo(sobol_examples)
.
Note that these procedures call fast C code. The following restrictions apply:
n
,d
must be .
n
must be and
d
must be .
if randomize = "none"
or randomize =
"digital.shift"
, n
must be and
d
must be .
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).
korobov()
and ghalton()
return an
-
matrix
; for an
-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)]))))) ## For sobol() examples, see demo(sobol_examples)
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)]))))) ## For sobol() examples, see demo(sobol_examples)
Functions for testing low-discrepancy sequences.
sum_of_squares(u) sobol_g(u, copula = copula::indepCopula(dim = ncol(u)), alpha = 1:ncol(u), ...) exceedance(x, q, p = 0.99, method = c("indicator", "individual.given.sum.exceeds", "sum.given.sum.exceeds"))
sum_of_squares(u) sobol_g(u, copula = copula::indepCopula(dim = ncol(u)), alpha = 1:ncol(u), ...) exceedance(x, q, p = 0.99, method = c("indicator", "individual.given.sum.exceeds", "sum.given.sum.exceeds"))
u |
|
copula |
|
alpha |
vector of parameters of Sobol's g test function. |
... |
additional arguments passed to the underlying
|
x |
|
q |
|
p |
If
|
method |
|
For examples see the demo man_test_functions
.
See ES_np(<matrix>)
from qrmtools for another test function.
sum_of_squares()
returns an -vector
(
numeric(n)
) with the rowwise computed scaled sum
of squares (theoretically integrating to 1).
sobol_g()
returns an -vector (
numeric(n)
)
with the rowwise computed Sobol' g functions.
exceedance()
's return value depends on method
:
returns indicators whether,
componentwise, x
exceeds the threshold determined by q
.
returns all rows of x
whose sum exceeds the threshold determined by q
.
returns the row sums of those
rows of x
whose sum exceeds the threshold determined by q
.
Marius Hofert and Christiane Lemieux
Radovic, I., Sobol', I. M. and Tichy, R. F. (1996). Quasi-Monte Carlo methods for numerical integration: Comparison of different low discrepancy sequences. Monte Carlo Methods and Applications 2(1), 1–14.
Faure, H., Lemieux, C. (2009). Generalized Halton Sequences in 2008: A Comparative Study. ACM-TOMACS 19(4), Article 15.
Owen, A. B. (2003). The dimension distribution and quadrature test functions. Stat. Sinica 13, 1-–17.
Sobol', I. M. and Asotsky, D. I. (2003). One more experiment on estimating high-dimensional integrals by quasi-Monte Carlo methods. Math. Comput. Simul. 62, 255–-263.
## Generate some (here: copula, pseudo-random) data library(copula) set.seed(271) cop <- claytonCopula(iTau(claytonCopula(), tau = 0.5)) # Clayton copula U <- rCopula(1000, copula = cop) ## Compute sum of squares test function mean(sum_of_squares(U)) # estimate of E(3(sum_{j=1}^d U_j^2)/d) ## Compute the Sobol' g test function if(packageVersion("copula") >= "0.999-20") mean(sobol_g(U)) # estimate of E(<Sobol's g function>) ## Compute an exceedance probability X <- qnorm(U) mean(exceedance(X, q = qnorm(0.99))) # fixed threshold q mean(exceedance(X, p = 0.99)) # empirically estimated marginal p-quantiles as thresholds ## Compute 99% expected shortfall for the sum mean(exceedance(X, p = 0.99, method = "sum.given.sum.exceeds")) ## Or use ES_np(X, level = 0.99) from 'qrmtools'
## Generate some (here: copula, pseudo-random) data library(copula) set.seed(271) cop <- claytonCopula(iTau(claytonCopula(), tau = 0.5)) # Clayton copula U <- rCopula(1000, copula = cop) ## Compute sum of squares test function mean(sum_of_squares(U)) # estimate of E(3(sum_{j=1}^d U_j^2)/d) ## Compute the Sobol' g test function if(packageVersion("copula") >= "0.999-20") mean(sobol_g(U)) # estimate of E(<Sobol's g function>) ## Compute an exceedance probability X <- qnorm(U) mean(exceedance(X, q = qnorm(0.99))) # fixed threshold q mean(exceedance(X, p = 0.99)) # empirically estimated marginal p-quantiles as thresholds ## Compute 99% expected shortfall for the sum mean(exceedance(X, p = 0.99, method = "sum.given.sum.exceeds")) ## Or use ES_np(X, level = 0.99) from 'qrmtools'
Converting higher-dimensional matrices of quasi-random numbers to arrays of specific formats.
to_array(x, f, format = c("(n*f,d)", "(n,f,d)"))
to_array(x, f, format = c("(n*f,d)", "(n,f,d)"))
x |
( |
f |
factor |
format |
|
to_array()
is helpful for converting quasi-random numbers
to time series paths.
(n * f, d)-matrix
or (n, f, d)-array
depending on the chosen format
.
Marius Hofert
korobov()
, ghalton()
, sobol()
.
## Basic call N <- 4 # replications n <- 3 # time steps d <- 2 # dimension set.seed(271) # note: respected for the choice of 'randomize' x <- sobol(N, d = n * d, randomize = "digital.shift") # higher-dim. Sobol' stopifnot(dim(to_array(x, f = n)) == c(N * n, d)) # conversion and check stopifnot(dim(to_array(x, f = n, format = "(n,f,d)")) == c(N, n, d)) ## See how the conversion is done (x <- matrix(1:(N * n * d), nrow = N, byrow = TRUE)) to_array(x, f = n) # => (n * d)-column x was blocked in n groups of size d each
## Basic call N <- 4 # replications n <- 3 # time steps d <- 2 # dimension set.seed(271) # note: respected for the choice of 'randomize' x <- sobol(N, d = n * d, randomize = "digital.shift") # higher-dim. Sobol' stopifnot(dim(to_array(x, f = n)) == c(N * n, d)) # conversion and check stopifnot(dim(to_array(x, f = n, format = "(n,f,d)")) == c(N, n, d)) ## See how the conversion is done (x <- matrix(1:(N * n * d), nrow = N, byrow = TRUE)) to_array(x, f = n) # => (n * d)-column x was blocked in n groups of size d each