mvTProbQuasiMonteCarlo {mvPot} | R Documentation |
Multivariate t distribution function
Description
Estimate the multivariate t distribution function with quasi-Monte Carlo method.
Usage
mvTProbQuasiMonteCarlo(p, upperBound, cov, nu, genVec, ...)
Arguments
p |
Number of samples used for quasi-Monte Carlo estimation. Must be a prime number. |
upperBound |
Vector of probabilities, i.e., the upper bound of the integral. |
cov |
Covariance matrix of the multivariate normal distribution. Must be positive semi-definite. WARNING: for performance in high-dimensions, no check is done to ensure positive-definiteness of the covariance matrix. It is the user responsibility to ensure that this property is verified. |
nu |
Degrees of freedom of the t distribution. |
genVec |
Generating vector for the quasi-Monte Carlo procedure. Can be computed using |
... |
Additional arguments passed to Cpp routine. |
Details
The function uses a quasi-Monte Carlo procedure based on randomly shifted lattice rules to estimate the distribution function a multivariate normal distribution as described in Genz and Bretz (2009) on page 50.
For compatibility reasons, the function handles the univariate case, which is passed on to pt
.
Value
A named vector with components estimate estimate
of the distribution function
along error
, 3 times the empirical Monte Carlo standard error over the nrep
replications.
Author(s)
Raphael de Fondeville
Source
Adapted from the QSILATMVTV
Matlab routine by A. Genz (2013)
References
Genz, A. and Bretz, F. (2009). Computations of Multivariate Normal and t Probabilities, volume 105. Springer: Dordrecht.
Examples
#Define locations
loc <- expand.grid(1:4, 1:4)
ref <- sample.int(16, 1)
#Define degrees of freedom
nu <- 3
#Compute variogram matrix
variogramMatrix <- ((sqrt((outer(loc[,1],loc[,1],"-"))^2 +
(outer(loc[,2],loc[,2],"-"))^2)) / 2)^(1.5)
#Define an upper bound
upperBound <- variogramMatrix[-ref,ref]
#Compute covariance matrix
cov <- (variogramMatrix[-ref,ref]%*%t(matrix(1, (nrow(loc) - 1), 1)) +
t(variogramMatrix[ref,-ref]%*%t(matrix(1, (nrow(loc) - 1), 1))) -
variogramMatrix[-ref,-ref])
#Compute generating vector
p <- 499
latticeRule <- genVecQMC(p, (nrow(loc) - 1))
#Estimate the multivariate distribution function
mvTProbQuasiMonteCarlo(latticeRule$primeP, upperBound, cov, nu, latticeRule$genVec)