mvtNormQuasiMonteCarlo {mvPot} | R Documentation |
Multivariate normal distribution function
Description
Estimate the multivariate distribution function with quasi-Monte Carlo method.
Usage
mvtNormQuasiMonteCarlo(p, upperBound, cov, 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 performed on the matrix. It is the user responsibility to ensure that the matrix is positive semi-definite. |
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.
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.
Source
Adapted from the QSILATMVNV
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)
#Compute variogram matrix
variogramMatrix <- ((sqrt((outer(loc[,1],loc[,1],"-"))^2 +
(outer(loc[,2],loc[,2],"-"))^2)) / 2)^(1.5)
#Define an upper boud
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
mvtNormQuasiMonteCarlo(latticeRule$primeP, upperBound, cov, latticeRule$genVec)