[R] Multivariate normal density in C for R

[Dimitris Rizopoulos]

On 6/26/2011 5:53 PM, zerfetzen wrote:
> IIRC, package mvtnorm will allow an X matrix, but requires mu to be a vector,
> so although it's close, it won't do it all...but all suggestions are well
> received.
>
> Dimitrius, you don't happen to have the multivariate t form of that
> function, do you?

Well, it's relatively easy to adjust it, e.g.,

dmvt <- function (x, mu, Sigma, df, log = FALSE) {
     if (!is.matrix(x))
         x <- rbind(x)
     p <- nrow(Sigma)
     ed <- eigen(Sigma, symmetric = TRUE)
     ev <- ed$values
     if (!all(ev >= -1e-06 * abs(ev[1])))
         stop("'Sigma' is not positive definite")
     ss <- if (!is.matrix(mu)) {
         x - rep(mu, each = nrow(x))
     } else {
         x - mu
     }
     inv.Sigma <- ed$vectors %*% (t(ed$vectors)/ev)
     quad <- rowSums((ss %*% inv.Sigma) * ss)/df
     fact <- lgamma((df + p)/2) - lgamma(df/2) -
         0.5 * (p * (log(pi) + log(df)) + sum(log(ev)))
     if (log)
         fact - 0.5 * (df + p) * log(1 + quad)
     else
         exp(fact) * ((1 + quad)^(-(df + p)/2))
}


Best,
Dimitris


> --
> View this message in context: http://r.789695.n4.nabble.com/Multivariate-normal-density-in-C-for-R-tp3624602p3626127.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at 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.
>

-- 
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014
Web: http://www.erasmusmc.nl/biostatistiek/