Skip to content
Prev 140109 / 398506 Next

Vectorization Problem

"Sergey Goriatchev" <sergeyg at gmail.com> wrote in
news:7cb007bd0803220542h66cfcd9awfd0a7a054d0fdaf8 at mail.gmail.com:
snipped your code that you didn't like

When Yan built his copula package, he called the dmvnorm function from 
Leisch's mvtnorm package:

dnormalCopula <- function(copula, u) {
  dim <- copula at dimension
  sigma <- getSigma(copula)
  if (is.vector(u)) u <- matrix(u, ncol = dim)
  x <- qnorm(u)
  val <- dmvnorm(x, sigma = sigma) / apply(x, 1, function(v) prod(dnorm
(v)))
  val[apply(u, 1, function(v) any(v <= 0))] <- 0
  val[apply(u, 1, function(v) any(v >= 1))] <- 0
  val
}

If the mvtnorm package is installed, one looks at the dmvnorm function 
simply by typing:

dmvnorm

I did not see any for-loops. After error checking, Leisch's code is:
--------
distval <- mahalanobis(x, center = mean, cov = sigma)
logdet <- sum(log(eigen(sigma, symmetric = TRUE, 
                         only.values = TRUE)$values))
logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
if (log) 
        return(logretval)
exp(logretval)
---------