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
rho <- 0.5 #corr
R <- rbind(c(1,rho),c(rho,1)) #vcov matrix
id <- diag(2)
gauss <- matrix(0,99,99)
u <- seq(0.01, 0.99, by=0.01)
for(i in 1:99){
for(j in 1:99){
xx <- t(t(c(qnorm(u[i]), qnorm(u[j]))))
gauss[i,j] <- 1/sqrt(1-rho^2)*exp(-0.5*t(xx)%*%(solve(R)-id)%*%xx)
}
}
persp(u, u, gauss, xlab="x1", ylab="x2", zlab="Density
of Gauss Copula",
theta=320, phi=20, col="palegreen",
ticktype="detailed")
--
Sergei
"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 Copulasnipped 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