Is there any vectorized equivalent of pmvnorm?

Question

I am trying to optimize my code which takes as argument four vectors X1, X2, X3, X4 and a correlation matrix R. First, for each observation I compute F(X3,X4)+F(X1,X2)-F(X2,X3)-F(X1,X4) where F corresponds to the function below:

F <- function(u,v,R){
pmvnorm(upper=c(qnorm(u,0,1), qnorm(v,0,1)),mean=c(0,0),lower=c(-Inf, -Inf),corr = R, sigma=NULL, algorithm = mvtnorm::GenzBretz(), keepAttr=FALSE))
}

Then I add up the (thresholded by a certain delta) logarithms of all the observations and return the opposite of the sum.

I would like my code to be faster and using mapply does not provide me with a significant difference as compared to a for loop.

Currently, I am using the code below:

L_n <- function(R, X1, X2, X3, X4){
  delta = 10**-9
  mysummands <- mapply(F, X3, X4, MoreArgs=list(R=R)) + mapply(F, X1, X2, MoreArgs=list(R=R)) - mapply(F, X2, X3, MoreArgs=list(R=R)) - mapply(F, X1, X4, MoreArgs=list(R=R))
  L <- sum(log(mapply(max,mysummands,MoreArgs=list(delta))))
  return(-L)
}

If I try using F as a vectorized function by specifying vectors for u and v, I get an error that the upper bound and the diagonal of the correlation matrix do not have the same size. (The correlation matrix R remains the same for all observations).

Does there exist a vectorized equivalent to pmvnorm ?

Answer 1

pmnorm in the mnormt package is vectorized:

library(mnormt)

mu <- c(1,12,2)
Sigma <- matrix(c(1,2,0,2,5,0.5,0,0.5,3), 3, 3)
pmnorm(c(2,11,3), mu, Sigma)
# [1] 0.2505509
pmnorm(c(3,10,2), mu, Sigma)
# [1] 0.1065538
pmnorm(
  rbind(
    c(2,11,3),
    c(3,10,2)
  ), mu, Sigma)
# [1] 0.2505509 0.1065538