dnvmix {nvmix} | R Documentation |
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 McNeil et al. (2015, Chapter 6). Supported are the following types of specification (see also the examples below):
|
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. |
factor |
(d, d) lower triangular |
control |
|
log |
|
verbose |
|
... |
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.
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.
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).
Erik Hintz, Marius Hofert and Christiane Lemieux.
McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.
### 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))