pnvmix {nvmix}R Documentation

Distribution Function of Multivariate Normal Variance Mixtures

Description

Evaluating multivariate normal variance mixture distribution functions (including Student t and normal distributions).

Usage


                   
pnvmix(upper, lower = matrix(-Inf, nrow = n, ncol = d), qmix,
      loc = rep(0, d), scale = diag(d), standardized = FALSE,
      control = list(), verbose = TRUE, ...)

pStudent(upper, lower = rep(-Inf, d), df, loc = rep(0, d), scale = diag(d),
         standardized = FALSE, control = list(), verbose = TRUE)
pNorm(upper, lower = rep(-Inf, d), loc = rep(0, d), scale = diag(d),
      standardized = FALSE, control = list(), verbose = TRUE)

Arguments

upper

(n, d)-matrix of upper integration limits; each row represents a d-vector of upper integration limits.

lower

(n, d)-matrix of lower integration limits (componentwise less than or equal to upper); each row represents a d-vector of lower integration limits.

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

standardized

logical indicating whether scale is assumed to be a correlation matrix.

control

list specifying algorithm specific parameters; see details below.

verbose

logical indicating whether a warning is given if the required precision abstol has not been reached; can also be an integer in which case 0 is FALSE, 1 is TRUE and 2 stands for producing a more verbose warning (for each set of provided integration bounds).

...

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

Details

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 increment) in each iteration until the pre-specified error tolerance abstol is reached. The attribute "numiter" gives the number of such iterations needed.

Algorithm specific parameters (such as the above mentioned 'increment' or 'abstol') 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.

mean.sqrt.mix

expectation of the square root √(W) of the mixing variable W. If NULL, it will be estimated via QMC; this is only needed for determining the reordering of the integration bounds, so a rather crude approximation is fine.

abstol

non-negative numeric providing the absolute precision required, defaults to 1e-3. If abstol = 0, the algorithm will typically run until the total number of function evaluations exceeds fun.eval[2]. If n > 1 (so x has more than one row), the algorithm runs until the precision requirement is reached for all n density estimates.

precond

logical indicating whether preconditioning is applied, that is, reordering of the integration variables. If TRUE, integration limits as well as scale are internally re-ordered in a way such that the overall variance of the integrand is usually smaller than with the original ordering; this usually leads smaller run-times.

increment

character string indicating how the sample size should be increased in each iteration. Available are:

"doubling":

next iteration has as many sample points as all the previous iterations combined.

"num.init":

all iterations use an additional fun.eval[1]-many points.

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, precond, 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 absolute error tolerance abstol cannot be achieved with fun.eval[2] function evaluations, an additional warning is thrown.

pStudent() and pNorm() are wrappers of pnvmix(, qmix = "inverse.gamma", df = df) and pnvmix(, qmix = "constant"), respectively. In the univariate case, the functions pt() and pnorm() are used.

Value

pnvmix(), pStudent() and pNorm() return a numeric n-vector with the computed probabilities and corresponding attributes "error" (error estimates of the RQMC estimator) and "numiter" (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.

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.

See Also

dnvmix(), rnvmix()

Examples

### 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 <- 0.5 # note that this is *non-integer*
set.seed(1)
pt1 <- pnvmix(b, lower = a, qmix = "inverse.gamma", df = df, scale = P)

## 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)
mean.sqrt.mix <- sqrt(df) * gamma(df/2) / (sqrt(2) * gamma((df+1) / 2))
set.seed(1)
pt2 <- pnvmix(b, lower = a, qmix = qW, 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(1)
pt3 <- pnvmix(b, lower = a, qmix = qW, df = df, scale = P)
stopifnot(all.equal(pt3, pt1, tol = 7e-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))

### 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 = 5e-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(1)
p.St.dfInf <- pStudent(b, df = Inf, scale = P)
set.seed(1)
p.Norm <- pNorm(b, scale = P)
stopifnot(all.equal(p.St.dfInf, p.Norm, check.attributes = FALSE))

[Package nvmix version 0.0-2 Index]