Re: [R] Multivariate normal density in C for R
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@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/ __ 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.
Re: [R] Multivariate normal density in C for R
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? -- 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@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.
Re: [R] Multivariate normal density in C for R
zerfetzen wrote: > > Does anyone know of a package that uses C code to calculate a multivariate > normal density? > > My goal is to find a faster way to calculate MVN densities and avoid R > loops or apply functions, such as when X and mu are N x K matrices, as > opposed to vectors, and in this particular case, speed really matters. I > would like to be able to use .C or .Call to pass X, mu, Sigma, and N to a > C program and have it return a vector of log densities to R. > > I'm new to putting C in R, but am sure I'll figure it out. Thanks for any > suggestions. > Just a suggestion: package mvtnorm Berend -- View this message in context: http://r.789695.n4.nabble.com/Multivariate-normal-density-in-C-for-R-tp3624602p3625954.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.
Re: [R] Multivariate normal density in C for R
Dimitris, Thanks for the great code. When the number of rows of X and mu are large, it is probably faster due to R's vectorization. Thanks again. -- View this message in context: http://r.789695.n4.nabble.com/Multivariate-normal-density-in-C-for-R-tp3624602p3625857.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.
Re: [R] Multivariate normal density in C for R
I use the following function which does not uses loops and seems to be pretty fast: dmvnorm <- 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 <- 0.5 * rowSums((ss %*% inv.Sigma) * ss) fact <- -0.5 * (p * log(2 * pi) + sum(log(ev))) if (log) as.vector(fact - quad) else as.vector(exp(fact - quad)) } I hope it helps. Best, Dimitris On 6/25/2011 3:58 PM, zerfetzen wrote: Does anyone know of a package that uses C code to calculate a multivariate normal density? My goal is to find a faster way to calculate MVN densities and avoid R loops or apply functions, such as when X and mu are N x K matrices, as opposed to vectors, and in this particular case, speed really matters. I would like to be able to use .C or .Call to pass X, mu, Sigma, and N to a C program and have it return a vector of log densities to R. I'm new to putting C in R, but am sure I'll figure it out. Thanks for any suggestions. -- View this message in context: http://r.789695.n4.nabble.com/Multivariate-normal-density-in-C-for-R-tp3624602p3624602.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. -- 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/ __ 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.