Density of Multivariate Normal Variance Mixtures
Evaluating multivariate normal variance mixture densities (including Student t and normal densities).
dnvmix(x, qmix, loc = rep(0, d), scale = diag(d),
factor = NULL, control = list(), log = FALSE, verbose = TRUE,...)
dStudent(x, df, loc = rep(0, d), scale = diag(d), factor = NULL,
log = FALSE, verbose = TRUE, ...)
dNorm(x, loc = rep(0, d), scale = diag(d), factor = NULL,
log = FALSE, verbose = TRUE, ...)x |
(n, d)- |
qmix |
specification of the mixing variable W; see
|
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); 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. Needs to be full rank for the density to exist. |
factor |
(d, d) lower triangular |
control |
|
log |
|
verbose |
|
... |
additional arguments (for example, parameters)
passed to the underlying mixing distribution when |
For the density to exist, scale must be full rank.
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().
dStudent() and dNorm() are wrappers of
dnvmix(, qmix = "inverse.gamma", df = df) and
dnvmix(, qmix = "constant"), respectively.
In these cases, dnvmix() uses the analytical formulas for the
density of a multivariate Student t and normal distribution,
respectively.
Internally, an adaptive randomized Quasi-Monte Carlo (RQMC) approach
is used to estimate the log-density. It is an iterative algorithm that
evaluates the integrand at a randomized Sobol' point-set (default) in
each iteration until the pre-specified error tolerance
control$dnvmix.reltol in the control argument
is reached for the log-density. The attribute
"numiter" gives the worst case number of such iterations needed
(over all rows of x). Note that this function calls underlying
C code.
Algorithm specific parameters (such as above mentioned control$dnvmix.reltol)
can be passed as a list via the control argument,
see get_set_param() for details and defaults.
If the error tolerance cannot be achieved within control$max.iter.rqmc
iterations and fun.eval[2] function evaluations, an additional
warning is thrown if verbose=TRUE.
dnvmix(), dStudent() and dNorm() return a
numeric n-vector with the computed (log-)density
values and attributes "abs. error" and "rel. error"
(containing the absolute and relative error
estimates of the of the (log-)density) and "numiter"
(containing the number of iterations).
Erik Hintz, Marius Hofert and Christiane Lemieux.
Hintz, E., Hofert, M. and Lemieux, C. (2020), Normal variance mixtures: Distribution, density and parameter estimation. https://arxiv.org/abs/1911.03017.
McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.
pnvmix(), rnvmix(), fitnvmix(),
get_set_param().
### Examples for dnvmix() ######################################################
## Generate a random correlation matrix in three dimensions
d <- 3
set.seed(271)
A <- matrix(runif(d * d), ncol = d)
P <- cov2cor(A %*% t(A))
## Evaluate a t_{3.5} density
df <- 3.5
x <- matrix(1:12/12, ncol = d) # evaluation points
dt1 <- dnvmix(x, qmix = "inverse.gamma", df = df, scale = P)
stopifnot(all.equal(dt1, c(0.013266542, 0.011967156, 0.010760575, 0.009648682),
tol = 1e-7, check.attributes = FALSE))
## Here is a version providing the quantile function of the mixing distribution
qW <- function(u, df) 1 / qgamma(1-u, shape = df/2, rate = df/2)
dt2 <- dnvmix(x, qmix = qW, df = df, scale = P)
## Compare
stopifnot(all.equal(dt1, dt2, tol = 5e-4, check.attributes = FALSE))
## Evaluate a normal density
dn <- dnvmix(x, qmix = "constant", scale = P)
stopifnot(all.equal(dn, c(0.013083858, 0.011141923, 0.009389987, 0.007831596),
tol = 1e-7, check.attributes = FALSE))
## Case with missing data
x. <- x
x.[3,2] <- NA
x.[4,3] <- NA
dt <- dnvmix(x., qmix = "inverse.gamma", df = df, scale = P)
stopifnot(is.na(dt) == rep(c(FALSE, TRUE), each = 2))
## Univariate case
x.. <- cbind(1:10/10) # (n = 10, 1)-matrix; vectors are considered rows in dnvmix()
dt1 <- dnvmix(x.., qmix = "inverse.gamma", df = df, factor = 1)
dt2 <- dt(as.vector(x..), df = df)
stopifnot(all.equal(dt1, dt2, check.attributes = FALSE))
### Examples for dStudent() and dNorm() ########################################
## Evaluate a t_{3.5} density
dt <- dStudent(x, df = df, scale = P)
stopifnot(all.equal(dt, c(0.013266542, 0.011967156, 0.010760575, 0.009648682),
tol = 1e-7, check.attributes = FALSE))
## Evaluate a normal density
x <- x[1,] # use just the first point this time
dn <- dNorm(x, scale = P)
stopifnot(all.equal(dn, 0.013083858, tol = 1e-7, check.attributes = FALSE))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.