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)
}
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