NLLmvnorm <- function(data, mu, sigma) {
  a <- qr(sigma)
  logdet <- sum(log(abs(diag(a$qr))))
  sigma.inv <- fast.ginv(sigma) #fast.ginv gives inverse of a matrix
  -0.5*logdet - 0.5*t(data-mu) %*% sigma.inv %*% (data-mu) - (length(data) / 2) * log(2 * pi) #MVN log likelihood