[R] Faster Multivariate Normal

[Duncan Murdoch]

On 07/06/2016 9:06 AM, Doran, Harold wrote:
> I am computing a complex multidimensional integral (4 dimensions) using Gauss-legendre and am looking to optimize one piece of code that is currently very expensive. The integral is evaluated for K individuals in my data and once this piece is computed it can be stored and recycled over all K individuals. Once this piece is stored, the integral can be computed in about 1.3 minutes using R for each individual.
>
> Because the integral has multiple dimensions, the number of nodes grows exponentially such that I require q^4 total nodes and experimentation is showing I need no fewer than 40 per dimension. At each of these, I need to compute the density of the multivariate normal at all q^4 nodes and store it. This is very expensive and I wonder if there is a way to improve the computational time to be faster?
>
> Below is just a reproducible toy example (not using legendre nodes) to illustrate the issue. The final line is what I am wondering about to see if it can be improved in terms of computational time.

I'd vectorize rather than using sapply (which is really a long loop). 
Try to put all the values into rows of a single matrix and just make one 
call to dmvnorm.

Duncan Murdoch

> Thank you
> Harold
>
> library(mvtnorm)
> ### Create parameters for MVN
> mu <- c(0,0,0,0)
> cov <- matrix(.2, ncol= 4,nrow=4)
> diag(cov) <- 1
> sigma <- as.matrix(cov)
>
> ### Create nodes and expand to 4 dimensions for quadrature
> dm  <- length(mu)
> gh  <- seq(-4, 4, length.out = 10)
> n <- length(gh)
> idx <- as.matrix(expand.grid(rep(list(1:n),dm)))
>
> ### Compute density, this is expensive
> adjFactor <- sapply(1:nrow(idx), function(i) dmvnorm(gh[idx[i,]], mean = mu, sigma = sigma))
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>