R/mvnpdf.R

In mavalosf/formationED:

#' mvnpdf multivariate normal density function
#'
#' multivariate normal density function of a p*n matrix x of mean vector mean and covariance matrix varcovM
#'
#' @param x p*n matrix of values
#' @param mean p vector of means
#' @param varcovM p*p symetric matrix of covariances
#' @param Log logic value default \code{TRUE} for logarithmic scale
#'
#' @return a list containing the input matrix x and y multivariate normal density values of x
#' @export
#'
#' @examples
#' mvnpdf(x=matrix(1.96),Log = F)
#' dnorm(1.96)
#' mvnpdf(x=matrix(1.96,ncol = 2),Log = F)
#'
mvnpdf <- function(x, mean =  rep(0, nrow(x)), varcovM = diag(nrow(x)), Log = TRUE) {

  n <- ncol(x)
  p <- nrow(x)

  x0 <- x - mean
  Rinv <- solve(varcovM)
  LogDetvarcovM <- log(det(varcovM))

  y<-apply(x0, 2, function(x){t(x) %*% Rinv %*% x})
  y<-p/2 * log(2*pi) - 0.5 * LogDetvarcovM -  0.5 * y

  if (!Log) {y <- exp(y)}

  solution<-list(x=x,y=y)
  return(solution)
}

# mvnpdf <- function(x, mean =  rep(0, nrow(x)),
#                    varcovM = diag(nrow(x)), Log = TRUE) {
#   n <- ncol(x)
#   p <- nrow(x)
#   x0 <- x - mean
#   Rinv <- somvnpdf(lve(varcovM)
#   LogDetvarcovM <- log(det(varcovM))
#
#   y <- NULL
#   for (j in 1:n) {
#     yj <- - p/2 * log(2*pi) - 0.5 * LogDetvarcovM -
#       0.5 * t(x0[, j]) %*% Rinv %*% x0[, j]
#     y <- c(y, yj)
#   }
#
#   if (!Log) {
#     y <- exp(y)
#   }
#
#   return(y)
# }