that's the bvnpdf() code: bvnpdf <- function(xy, mu, Siginv, detSig){ f<-numeric() x <- xy[[1]] if(xy[[1]]>xy[[2]] & (1==2)) { f <- 0} else { v <- rbind(xy[1], xy[2]) e <- t(v-mu) %*% Siginv %*% (v-mu) f <- as.numeric(exp(-e/2)/(2*pi)/sqrt(detSig)) } f }