From: David Winsemius <dwinsemius_at_comcast.net>

Date: Sat, 22 Mar 2008 16:59:57 +0000 (UTC)

distval <- mahalanobis(x, center = mean, cov = sigma) logdet <- sum(log(eigen(sigma, symmetric = TRUE,

Date: Sat, 22 Mar 2008 16:59:57 +0000 (UTC)

"Sergey Goriatchev" <sergeyg_at_gmail.com> wrote in news:7cb007bd0803220542h66cfcd9awfd0a7a054d0fdaf8_at_mail.gmail.com:

> I have the code for the bivariate Gaussian copula. It is written

*> with for-loops, it works, but I wonder if there is a way to
**> vectorize the function.
**> I don't see how outer() can be used in this case, but maybe one can
**> use mapply() or Vectorize() in some way? Could anyone help me,
**> please?
**>
**> ## Density of Gauss Copula
*

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)

-- David Winsemius ______________________________________________ R-help_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.Received on Sat 22 Mar 2008 - 18:42:00 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sat 22 Mar 2008 - 22:30:24 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*