Distribution Function of Multivariate Normal Variance Mixtures
Evaluating multivariate normal variance mixture distribution functions (including Student t and normal distributions).
pnvmix(upper, lower = matrix(-Inf, nrow = n, ncol = d), qmix, rmix,
loc = rep(0, d), scale = diag(d), standardized = FALSE,
control = list(), verbose = TRUE, ...)
pStudent(upper, lower = matrix(-Inf, nrow = n, ncol = d), df, loc = rep(0, d),
scale = diag(d), standardized = FALSE, control = list(), verbose = TRUE)
pNorm(upper, lower = matrix(-Inf, nrow = n, ncol = d), loc = rep(0, d),
scale = diag(d), standardized = FALSE, control = list(), verbose = TRUE)upper |
(n, d)- |
lower |
(n, d)- |
qmix, rmix |
specification of the mixing variable W via a quantile
function (
|
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. |
standardized |
|
control |
|
verbose |
|
... |
additional arguments (for example, parameters) passed
to the underlying mixing distribution when |
One should highlight that evaluating normal variance mixtures is a non-trivial tasks which, at the time of development of nvmix, was not available in R before, not even the special case of a multivariate Student t distribution for non-integer degrees of freedom, which frequently appears in applications in finance, insurance and risk management after estimating such distributions.
Note that the procedures call underlying C code. Currently, dimensions
d >= 16510 are not supported for the default method
sobol.
Internally, an iterative randomized Quasi-Monte Carlo (RQMC) approach
is used to estimate the probabilities. It is an iterative algorithm
that evaluates the integrand at a point-set (with size as specified by
control$increment in the control argument) in each
iteration until the pre-specified absolute error tolerance
control$pnvmix.abstol (or relative error tolerance
control$pnvmix.reltol which is used only when
control$pnvmix.abstol = NA) is reached. The attribute
"numiter" gives the number of such iterations needed.
Algorithm specific parameters (such as the above mentioned
control$pnvmix.abstol) can be passed as a list
via control, see get_set_param() for more
details. If specified error tolerances are not reached and
verbose = TRUE, a warning is thrown.
If provided scale
is singular, pnvmix() estimates the correct probability but
throws a warning if verbose = TRUE.
It is recommended to supply a quantile function via qmix, if available,
as in this case efficient RQMC methods are used to approximate the probability.
If rmix is provided, internally used is
plain MC integration, typically leading to slower convergence.
If both qmix and rmix are provided, the latter
is ignored.
pnvmix(), pStudent() and pNorm() return a
numeric n-vector with the computed probabilities
and corresponding attributes "abs. error" and rel. error
(error estimates of the RQMC estimator) and "numiter" (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.
Genz, A. and Bretz, F. (1999). Numerical computation of multivariate t-probabilities with application to power calculation of multiple contrasts. Journal of Statistical Computation and Simulation 63(4), 103–117.
Genz, A. and Bretz, F. (2002). Comparison of methods for the computation of multivariate t probabilities. Journal of Computational and Graphical Statistics 11(4), 950–971.
Genz, A. and Kwong, K. (2000). Numerical evaluation of singular multivariate normal distributions. Journal of Statistical Computation and Simulation 68(1), 1–21.
dnvmix(), rnvmix(), fitnvmix(),
pgnvmix(), get_set_param()
### Examples for pnvmix() ######################################################
## 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))
## Evaluate a t_{1/2} distribution function
a <- -3 * runif(d) * sqrt(d) # random lower limit
b <- 3 * runif(d) * sqrt(d) # random upper limit
df <- 1.5 # note that this is *non-integer*
set.seed(123)
pt1 <- pnvmix(b, lower = a, qmix = "inverse.gamma", df = df, scale = P)
## Here is a version providing the quantile function of the mixing distribution
qmix <- function(u, df) 1 / qgamma(1-u, shape = df/2, rate = df/2)
mean.sqrt.mix <- sqrt(df) * gamma(df/2) / (sqrt(2) * gamma((df+1) / 2))
set.seed(123)
pt2 <- pnvmix(b, lower = a, qmix = qmix, df = df, scale = P,
control = list(mean.sqrt.mix = mean.sqrt.mix))
## Compare
stopifnot(all.equal(pt1, pt2, tol = 7e-4, check.attributes = FALSE))
## mean.sqrt.mix will be approximated by QMC internally if not provided,
## so the results will differ slightly.
set.seed(123)
pt3 <- pnvmix(b, lower = a, qmix = qmix, df = df, scale = P)
stopifnot(all.equal(pt3, pt1, tol = 7e-4, check.attributes = FALSE))
## Here is a version providing a RNG for the mixing distribution
## Note the significantly larger number of iterations in the attribute 'numiter'
## compared to when 'qmix' was provided (=> plain MC versus RQMC)
set.seed(123)
pt4 <- pnvmix(b, lower = a,
rmix = function(n, df) 1/rgamma(n, shape = df/2, rate = df/2),
df = df, scale = P)
stopifnot(all.equal(pt4, pt1, tol = 8e-4, check.attributes = FALSE))
## Case with missing data and a matrix of lower and upper bounds
a. <- matrix(rep(a, each = 4), ncol = d)
b. <- matrix(rep(b, each = 4), ncol = d)
a.[2,1] <- NA
b.[3,2] <- NA
pt <- pnvmix(b., lower = a., qmix = "inverse.gamma", df = df, scale = P)
stopifnot(is.na(pt) == c(FALSE, TRUE, TRUE, FALSE))
## Case where upper = (Inf,..,Inf) and lower = (-Inf,...,-Inf)
stopifnot(all.equal(pnvmix(upper = rep(Inf, d), qmix = "constant"), 1,
check.attributes = FALSE))
## An example with singular scale:
A <- matrix( c(1, 0, 0, 0,
2, 1, 0, 0,
3, 0, 0, 0,
4, 1, 0, 1), ncol = 4, nrow = 4, byrow = TRUE)
scale <- A%*%t(A)
upper <- 2:5
pn <- pnvmix(upper, qmix = "constant", scale = scale) # multivariate normal
pt <- pnvmix(upper, qmix = "inverse.gamma", scale = scale, df = df) # multivariate t
stopifnot(all.equal(pn, 0.8581, tol = 1e-3, check.attributes = FALSE))
stopifnot(all.equal(pt, 0.7656, tol = 1e-3, check.attributes = FALSE))
## Evaluate a Exp(1)-mixture
## Specify the mixture distribution parameter
rate <- 1.9 # exponential rate parameter
## Method 1: Use R's qexp() function and provide a list as 'mix'
set.seed(42)
(p1 <- pnvmix(b, lower = a, qmix = list("exp", rate = rate), scale = P))
## Method 2: Define the quantile function manually (note that
## we do not specify rate in the quantile function here,
## but conveniently pass it via the ellipsis argument)
set.seed(42)
(p2 <- pnvmix(b, lower = a, qmix = function(u, lambda) -log(1-u)/lambda,
scale = P, lambda = rate))
## Check
stopifnot(all.equal(p1, p2))
### Examples for pStudent() and pNorm() ########################################
## Evaluate a t_{3.5} distribution function
set.seed(271)
pt <- pStudent(b, lower = a, df = 3.5, scale = P)
stopifnot(all.equal(pt, 0.6180, tol = 7e-5, check.attributes = FALSE))
## Evaluate a normal distribution function
set.seed(271)
pn <- pNorm(b, lower = a, scale = P)
stopifnot(all.equal(pn, 0.7001, tol = 1e-4, check.attributes = FALSE))
## pStudent deals correctly with df = Inf:
set.seed(123)
p.St.dfInf <- pStudent(b, df = Inf, scale = P)
set.seed(123)
p.Norm <- pNorm(b, scale = P)
stopifnot(all.equal(p.St.dfInf, p.Norm, check.attributes = FALSE))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.