rnvmix | R Documentation |
Generate vectors of random variates from multivariate normal variance mixtures (including Student t and normal distributions).
rnvmix(n, rmix, qmix, loc = rep(0, d), scale = diag(2), factor = NULL, method = c("PRNG", "sobol", "ghalton"), skip = 0, ...) rStudent(n, df, loc = rep(0, d), scale = diag(2), factor = NULL, method = c("PRNG", "sobol", "ghalton"), skip = 0) rNorm(n, loc = rep(0, d), scale = diag(2), factor = NULL, method = c("PRNG", "sobol", "ghalton"), skip = 0) rNorm_sumconstr(n, weights, s, method = c("PRNG", "sobol", "ghalton"), skip = 0)
n |
sample size n (positive integer). |
rmix |
specification of the mixing variable W, see McNeil
et al. (2015, Chapter 6) and Hintz et al. (2020),
via a random number generator.
This argument is ignored for
|
qmix |
specification of the mixing variable W via a quantile
function. This argument is required
for
|
df |
positive degress of freedom; can also be |
loc |
location vector of dimension d; this equals the mean vector of a random vector following the specified normal variance mixture distribution if and only if the latter exists. |
scale |
scale matrix (a covariance matrix entering the
distribution as a parameter) of dimension (d, d) (defaults to
d = 2);
this equals the covariance matrix of a random vector following
the specified normal variance mixture distribution divided by
the expecation of the mixing variable W if and only if the
former exists.
Note that |
factor |
(d, k)- |
method |
If |
skip |
|
weights |
d- |
s |
|
... |
additional arguments (for example, parameters) passed to
the underlying mixing distribution when |
Internally used is factor
, so scale
is not required
to be provided if factor
is given.
The default factorization used to obtain factor
is the Cholesky
decomposition via chol()
. To this end, scale
needs to have full rank.
Sampling from a singular normal variance mixture distribution can be
achieved by providing factor
.
The number of rows of factor
equals the dimension d of
the sample. Typically (but not necessarily), factor
is square.
rStudent()
and rNorm()
are wrappers of
rnvmix(, qmix = "inverse.gamma", df = df)
and
rnvmix(, qmix = "constant", df = df)
, respectively.
The function rNorm_sumconstr()
can be used to sample from the
multivariate standard normal distribution under a weighted sum constraint;
the implementation is based on Algorithm 1 in Vrins (2018). Let
Z = (Z_1,…,Z_d)~N_d(0, I_d). The function rNorm_sumconstr()
then samples from Z | w^T Z = s where w and s correspond
to the arguments weights
and s
. If supplied s
is a vector
of length n
, the i'th row of the returned matrix uses the constraint
w^T Z = s_i where s_i is the i'th element in s
.
rnvmix()
returns an (n, d)-matrix
containing n samples of the specified (via mix
)
d-dimensional multivariate normal variance mixture with
location vector loc
and scale matrix scale
(a covariance matrix).
rStudent()
returns samples from the d-dimensional
multivariate Student t distribution with location vector
loc
and scale matrix scale
.
rNorm()
returns samples from the d-dimensional
multivariate normal distribution with mean vector
loc
and covariance matrix scale
.
rNorm_sumconstr()
returns samples from the d-dimensional
multivariate normal distribution conditional on the weighted sum being
constrained to s
.
Erik Hintz, Marius Hofert and Christiane Lemieux
Hintz, E., Hofert, M. and Lemieux, C. (2021), Normal variance mixtures: Distribution, density and parameter estimation. Computational Statistics and Data Analysis 157C, 107175.
Hintz, E., Hofert, M. and Lemieux, C. (2022), Multivariate Normal Variance Mixtures in R: The R Package nvmix. Journal of Statistical Software, doi: 10.18637/jss.v102.i02.
McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.
Vrins, E. (2018) Sampling the Multivariate Standard Normal Distribution under a Weighted Sum Constraint. Risks 6(3), 64.
dnvmix()
, pnvmix()
### Examples for rnvmix() ###################################################### ## Generate a random correlation matrix in d dimensions d <- 3 set.seed(157) A <- matrix(runif(d * d), ncol = d) P <- cov2cor(A %*% t(A)) ## Draw random variates and compare df <- 3.5 n <- 1000 set.seed(271) X <- rnvmix(n, rmix = "inverse.gamma", df = df, scale = P) # with scale set.seed(271) X. <- rnvmix(n, rmix = "inverse.gamma", df = df, factor = t(chol(P))) # with factor stopifnot(all.equal(X, X.)) ## Checking df = Inf set.seed(271) X <- rnvmix(n, rmix = "constant", scale = P) # normal set.seed(271) X. <- rnvmix(n, rmix = "inverse.gamma", scale = P, df = Inf) # t_infinity stopifnot(all.equal(X, X.)) ## Univariate case (dimension = number of rows of 'factor' = 1 here) set.seed(271) X.1d <- rnvmix(n, rmix = "inverse.gamma", df = df, factor = 1/2) set.seed(271) X.1d. <- rnvmix(n, rmix = "inverse.gamma", df = df, factor = 1)/2 # manual scaling stopifnot(all.equal(X.1d, X.1d.)) ## Checking different ways of providing 'mix' ## 1) By providing a character string (and corresponding ellipsis arguments) set.seed(271) X.mix1 <- rnvmix(n, rmix = "inverse.gamma", df = df, scale = P) ## 2) By providing a list; the first element has to be an existing distribution ## with random number generator available with prefix "r" rinverse.gamma <- function(n, df) 1 / rgamma(n, shape = df/2, rate = df/2) set.seed(271) X.mix2 <- rnvmix(n, rmix = list("inverse.gamma", df = df), scale = P) ## 3) The same without extra arguments (need the extra list() here to ## distinguish from Case 1)) rinverseGamma <- function(n) 1 / rgamma(n, shape = df/2, rate = df/2) set.seed(271) X.mix3 <- rnvmix(n, rmix = list("inverseGamma"), scale = P) ## 4) By providing a quantile function ## Note: P(1/Y <= x) = P(Y >= 1/x) = 1-F_Y(1/x) = y <=> x = 1/F_Y^-(1-y) set.seed(271) X.mix4 <- rnvmix(n, qmix = function(p) 1/qgamma(1-p, shape = df/2, rate = df/2), scale = P) ## 5) By providing random variates set.seed(271) # if seed is set here, results are comparable to the above methods W <- rinverse.gamma(n, df = df) X.mix5 <- rnvmix(n, rmix = W, scale = P) ## Compare (note that X.mix4 is not 'all equal' with X.mix1 or the other samples) ## since rgamma() != qgamma(runif()) (or qgamma(1-runif())) stopifnot(all.equal(X.mix2, X.mix1), all.equal(X.mix3, X.mix1), all.equal(X.mix5, X.mix1)) ## For a singular normal variance mixture: ## Need to provide 'factor' A <- matrix( c(1, 0, 0, 1, 0, 1), ncol = 2, byrow = TRUE) stopifnot(all.equal(dim(rnvmix(n, rmix = "constant", factor = A)), c(n, 3))) stopifnot(all.equal(dim(rnvmix(n, rmix = "constant", factor = t(A))), c(n, 2))) ## Using 'skip'. Need to reset the seed everytime to get the same shifts in "sobol". ## Note that when using method = "sobol", we have to provide 'qmix' instead of 'rmix'. set.seed(271) X.skip0 <- rnvmix(n, qmix = "inverse.gamma", df = df, scale = P, method = "sobol") set.seed(271) X.skip1 <- rnvmix(n, qmix = "inverse.gamma", df = df, scale = P, method = "sobol", skip = n) set.seed(271) X.wo.skip <- rnvmix(2*n, qmix = "inverse.gamma", df = df, scale = P, method = "sobol") X.skip <- rbind(X.skip0, X.skip1) stopifnot(all.equal(X.wo.skip, X.skip)) ### Examples for rStudent() and rNorm() ######################################## ## Draw N(0, P) random variates by providing scale or factor and compare n <- 1000 set.seed(271) X.n <- rNorm(n, scale = P) # providing scale set.seed(271) X.n. <- rNorm(n, factor = t(chol(P))) # providing the factor stopifnot(all.equal(X.n, X.n.)) ## Univariate case (dimension = number of rows of 'factor' = 1 here) set.seed(271) X.n.1d <- rNorm(n, factor = 1/2) set.seed(271) X.n.1d. <- rNorm(n, factor = 1)/2 # manual scaling stopifnot(all.equal(X.n.1d, X.n.1d.)) ## Draw t_3.5 random variates by providing scale or factor and compare df <- 3.5 n <- 1000 set.seed(271) X.t <- rStudent(n, df = df, scale = P) # providing scale set.seed(271) X.t. <- rStudent(n, df = df, factor = t(chol(P))) # providing the factor stopifnot(all.equal(X.t, X.t.)) ## Univariate case (dimension = number of rows of 'factor' = 1 here) set.seed(271) X.t.1d <- rStudent(n, df = df, factor = 1/2) set.seed(271) X.t.1d. <- rStudent(n, df = df, factor = 1)/2 # manual scaling stopifnot(all.equal(X.t.1d, X.t.1d.)) ## Check df = Inf set.seed(271) X.t <- rStudent(n, df = Inf, scale = P) set.seed(271) X.n <- rNorm(n, scale = P) stopifnot(all.equal(X.t, X.n)) ### Examples for rNorm_sumconstr() ############################################# set.seed(271) weights <- c(1, 1) Z.constr <- rNorm_sumconstr(n, weights = c(1, 1), s = 2) stopifnot(all(rowSums(Z.constr ) == 2)) plot(Z.constr , xlab = expression(Z[1]), ylab = expression(Z[2]))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.