How do you code the likelihood of the normal distribution kernel in R manually?

Specifically, how do you code the product of the differences of x and mu, precision matrix, and transpose of the differences of x and mu. Is my code below correct? Thanks in advance.

(colSums(dat-mu_mat)%*%solve(sigma)%*%colSums(dat-mu_mat))

where mu_mat is the row vector of means repeated down n times.

Full code is below:

dat<-rmvnorm(100,mean=c(200,0.1),sigma=matrix(c(5,0,0,0.02),nrow=2))
n<-nrow(dat)
mu<-matrix(c(200,0.1),nrow=1)
mu_mat<-matrix(rep(c(200,0.1),100),nrow=100,ncol=2,byrow=TRUE)

loglik_mvn<-function(n,d,x,mu_mat,sigma){
  (-n*d/2)*log(2*pi)-(n/2)*det(sigma,log=TRUE)-0.5*(colSums(x-mu_mat)%*%solve(sigma)%*%colSums(x-mu_mat))

loglik_mvn(nrow(dat),ncol(dat),dat,mu_mat,sigma)
dmvnorm(dat,mean=mu,sigma=sigma,log=TRUE)
}

Solution

if you have dataset and a vector of means, and the sigma, then you could compute the exponent part as:

sum(diag(solve(sig, tcrossprod(t(x) - mu))))

Therefore you likelihood will look like:

log_like <- function(x, mu, sig){
  -sum(diag(solve(sig, tcrossprod(t(x) - mu))))/2 - #The part you are interested in
    nrow(x)*(ncol(x)*log(2*pi) + log(det(sig)))/2
}

x <- iris[-5]
mu <- colMeans(x)
sig <- var(x)

log_like(x, mu, sig)
#> [1] -379.9213

sum(mvtnorm::dmvnorm(x, mu, sig, TRUE)) # Produces same results as above
#> [1] -379.9213

^{Created on 2023-02-16 with reprex v2.0.2}

You could also use the following:

log_like2 <- function(x, mu, sig){
  u <- t(x) - mu #Notice mu is a vector and not a matrix
  -sum(u*solve(sig, u))/2 -  # the part you are interested in
   nrow(x)*(ncol(x)*log(2*pi) + log(det(sig)))/2
}

log_like2(x, mu, sig)
#> [1] -379.9213

putting all things together and using the recycling of R:

log_like3 <- function(x, mu, sig){
  u <- t(x) - mu
  -sum(u * solve(sig, u) + log(2*pi) + log(det(sig))/ncol(x))/2
}
log_like3(x, mu, sig)
[1] -379.9213