### five observations
N <- 5L
### dimension
J <- 4L
### lower and upper bounds, ie interval-censoring
lwr <- matrix(-runif(N * J), nrow = J)
upr <- matrix(runif(N * J), nrow = J)
### Cholesky factor
(C <- ltMatrices(runif(J * (J + 1) / 2), diag = TRUE))
### corresponding covariance matrix
(S <- as.array(Tcrossprod(C))[,,1])
### plain Monte-Carlo (Genz, 1992)
w <- NULL
M <- 25000
### quasi-Monte-Carlo (Genz & Bretz, 2002, but with different weights)
if (require("qrng")) w <- t(ghalton(M * N, J - 1))
### log-likelihood
lpmvnorm(lower = lwr, upper = upr, chol = C, w = w, M = M)
### compare with pmvnorm
exp(lpmvnorm(lower = lwr, upper = upr, chol = C, logLik = FALSE, w = w, M = M))
sapply(1:N, function(i) pmvnorm(lower = lwr[,i], upper = upr[,i], sigma = S))
### log-lik contributions and score matrix
slpmvnorm(lower = lwr, upper = upr, chol = C, w = w, M = M, logLik = TRUE)
Run the code above in your browser using DataLab