dnvmix {nvmix}R Documentation

Density of Multivariate Normal Variance Mixtures

Description

Evaluating multivariate normal variance mixture densities (including Student t and normal densities).

Usage

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, ...)

Arguments

x

(n, d)-matrix of evaluation points.

qmix

specification of the mixing variable W; see McNeil et al. (2015, Chapter 6). Supported are the following types of specification (see also the examples below):

character:

character string specifying a supported mixing distribution; currently available are "constant" (in which case W = 1 and thus the multivariate normal distribution with mean vector loc and covariance matrix scale results) and "inverse.gamma" (in which case W is inverse gamma distributed with shape and rate parameters df/2 and thus the multivariate Student t distribution with df degrees of freedom (required to be provided via the ellipsis argument) results).

list:

list of length at least one, where the first component is a character string specifying the base name of a distribution whose quantile function can be accessed via the prefix "q"; an example is "exp" for which "qexp" exists. If the list is of length larger than one, the remaining elements contain additional parameters of the distribution; for "exp", for example, this can be the parameter rate.

function:

function interpreted as the quantile function of the mixing variable W.

df

positive degress of freedom; can also be Inf in which case the distribution is interpreted as the multivariate normal distribution with mean vector loc and covariance matrix scale).

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.

factor

(d, d) lower triangular matrix such that factor %*% t(factor) equals scale; note that for performance reasons, this property is not tested. If not provided, factor is internally determined via t(chol()).

control

list specifying algorithm specific parameters; see details below.

log

logical indicating whether the logarithmic density is to be computed.

verbose

logical indicating whether a warning is given if the required precision abstol has not been reached.

...

additional arguments (for example, parameters) passed to the underlying mixing distribution when qmix is a character string or function.

Details

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.

The number of rows of factor equals the dimension d of the sample. Typically (but not necessarily), factor is square.

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 iterative randomized Quasi-Monte Carlo (RQMC) approach is used to estimate the 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 abstol is reached for both the density and 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 can be passed as a list via control. It can contain any of the following:

method

character string indicating the method to be used to compute the integral. Available are:

"sobol":

Sobol' sequence (default).

"ghalton":

generalized Halton sequence.

"PRNG":

plain Monte Carlo based on a pseudo-random number generator.

dnvmix.reltol

non-negative numeric providing the relative precision required, defaults to 0.025. If not provided, dnvmix.abstol will be used.

dnvmix.abstol

non-negative numeric providing the absolute precision required. Used only when dnvmix.reltol is NA.

CI.factor

multiplier of the Monte Carlo confidence interval bounds. The algorithm runs until CI.factor times the estimated standard error is less than abstol. If CI.factor = 3.3 (the default), one can expect the actual absolute error to be less than abstol in 99.9% of the cases.

fun.eval

numeric(2) providing the size of the first point set to be used to estimate the densities (typically a power of 2) and the maximal number of function evaluations. fun.eval defaults to c(2^7, 1e8).

max.iter.rqmc

numeric, providing the maximum number of iterations allowed in the RQMC approach; the default is 15.

B

number of randomizations for obtaining an error estimate in the randomized quasi-Monte Carlo (RQMC) approach; the default is 12.

Care should be taken when changing the algorithm-specific parameters, notably method, fun.eval[2] and B. Error estimates will not be reliable for too small B and the performance of the algorithm depends heavily on the (quasi-)Monte Carlo point-set used.

If the error tolerance dnvmix.reltol (or, if not supplied, dnvmix.reltol) cannot be achieved within max.iter.rqmc iterations and fun.eval[2] function evaluations, an additional warning is thrown.

Value

dnvmix(), dStudent() and dNorm() return a numeric n-vector with the computed (log-)density values and attributes "error" (containing the absolte error estimates of the of the (log-)density) and "numiter" (containing the number of iterations).

Author(s)

Erik Hintz, Marius Hofert and Christiane Lemieux.

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

See Also

pnvmix(), rnvmix()

Examples

### 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; note: vectors are taken as 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))

[Package nvmix version 0.0-2 Index]