mvTProbQuasiMonteCarlo | R Documentation |
Estimate the multivariate t distribution function with quasi-Monte Carlo method.
mvTProbQuasiMonteCarlo(p, upperBound, cov, nu, genVec, ...)
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. |
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
.
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.
Raphael de Fondeville
Adapted from the QSILATMVTV
Matlab routine by A. Genz (2013)
Genz, A. and Bretz, F. (2009). Computations of Multivariate Normal and t Probabilities, volume 105. Springer: Dordrecht.
#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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.