I also tried:
m1 <- length(x)-1 X1<- cbind(x[1:m1],x1[2:length(x)]) X2<- cbind(x[1:m1],x1[2:length(x)]) integral <- function(rho){ m1 <- length(x1)-1 integral <- apply(X2,1,function(y) apply(X1,1,function(x) pmvnorm(lower=c(y[1],x[1]), upper=c( y[2], x[2]), corr=matrix(c(1,rho,rho,1),ncol=2))) ) return(integral)} integral(0.5) But it's even slower. I thougt apply function might improve speed, maybe I'm using it wrong. In the current status it would take weeks to calculate, because I have to run it nearly 500.000 thousand times. -- View this message in context: http://r.789695.n4.nabble.com/Bivariate-normal-integral-tp4571018p4571362.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help@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.