pnvmix {nvmix} | R Documentation |
Evaluating multivariate normal variance mixture distribution functions (including Student t and normal distributions).
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)
upper |
(n, d)- |
lower |
(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. |
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 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.
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).
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.
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.
### 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))