mvtNormQuasiMonteCarlo {mvPot}

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

...

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.

References

Genz, A. and Bretz, F. (2009). Computations of Multivariate Normal and t Probabilities, volume 105. Springer: Dordrecht.

Genz, A. (2013). QSILATMVNV http://www.math.wsu.edu/faculty/genz/software/software.html

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)