# R/dmvnorm.R In csgillespie/abcpmcmc: Scalable Inference for Markov Processes with Intractable Likelihoods

```##################################################################
# Public functions
##################################################################
#' Multivariate Normal Density
#'
#' Calculates mutlivariate normal densities. This is different to the dmvnorm function
#' in the mvtnorm package in that it takes a matrix for both x and mean. It then calculates
#' a vector of densities according to dmvnorm(x[i,],mean[i,],sigma,log = FALSE).
#' To aid computation the mahalanobis distances are calculated in parallel using mclapply.
#' @importFrom parallel mclapply
#' @param x A matrix of values
#' @param mean A matrix of means
#' @param sigma A covariance matrix
#' @param log Boolean for whether we want log densities or not
#' @details My own implementation of the multivariate normal density function
#' for increased efficiency for application in this package
#' because there are so many repeated calls to densitymvnorm using a
#' given sigma matrix it made sense to have one that could take a matrix
#' of means as well as a matrix of x's and treat them as
#' paired, returning the density of x[1,] given mean[1,], x[2,]
#' given mean[2,] also stops repeated inversions of the matrix sigma, and
#' calculates the densities in parallel
#' @return A vector of densities
#' @export
densityMvNorm = function (x, mean, sigma, log = FALSE)
{
## takes a matrix of means rather than a single vector
if (missing(sigma)) sigma = diag(ncol(x))

if (NCOL(x) != NCOL(sigma)) {
stop("x and sigma have non-conforming size")
}
if (!isSymmetric(sigma, tol = sqrt(.Machine\$double.eps),
check.attributes = FALSE)) {
stop("sigma must be a symmetric matrix")
}
if (NCOL(mean) != NROW(sigma)) {
stop("mean and sigma have non-conforming size")
}
## invert matrix before hand so only do this once
prec = solve(sigma)
means = lapply(1:dim(mean)[1],function(i){mean[i,]})
distval = do.call(rbind,
mclapply(means,
mahalanobis,x = x, cov = prec,inverted = TRUE))
logdet = sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)\$values))
logretval = -(ncol(x) * log(2 * pi) + logdet + distval)/2

if (log)
return(logretval)
else
return(exp(logretval))
}
```
csgillespie/abcpmcmc documentation built on May 14, 2019, 12:11 p.m.