Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

The formula for the multivariate normal density function in the standard form contains $\Sigma^{-1}$ and the determinant of $\Sigma$, which are not very computationally friendly. Is it possible to calculate the density in some other way, that is more suitable for computer implementation?

EDIT: To be precise I am after the log of the density function, if it helps.

share|improve this question
    
why do you say $\Sigma$ is not computer friendly? –  mathemagician Jan 17 '13 at 13:13
    
@mathemagician, calculating determinants and inverse matrices is not computer friendly in general, because they are slow and numerically unstable (that's what I have been told at least!) –  Grzenio Jan 17 '13 at 13:15
    
It is not that bad if you use good software packages. It could all go wrong tho if the matrix has both very large and very small numbers, as a process called floating point numbers is used. I don't know any further about this, but I use R and it works fine, even with very large matrices. –  mathemagician Jan 17 '13 at 13:25
add comment

1 Answer

Perhaps this can help you:

http://gallery.rcpp.org/articles/dmvnorm_arma/

They use an alternative implementation of the function dmvnorm() from package mvtnorm() in R. It's a quite good time reduction, take a look at this example I just ran:

> x=matrix(0,1,40)
> m=rep(0,40)
> s=diag(40)
> dmvnorm(x,m,s)
[1] 1.087433e-16
> dMvn(x,m,s)
[1] 1.087433e-16
> benchmark(dmvnorm(x,m,s), dMvn(x,m,s), replications=10000)[,1:4]
              test replications elapsed relative
2    dMvn(x, m, s)        10000    1.35    1.000
1 dmvnorm(x, m, s)        10000    5.07    3.756

I have a C implementation of dmvnorm that is taking more than 95% of the total running time. Now, I will implement this dMvn() and I hope that the bottleneck can tone down. The function is:

dMvn <- function(X,mu,Sigma) {
    k <- ncol(X)
    rooti <- backsolve(chol(Sigma),diag(k))
    quads <- colSums((crossprod(rooti,(t(X)-mu)))^2)
    return(exp(-(k/2)*log(2*pi) + sum(log(diag(rooti))) - .5*quads))
}
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.