calculate density of multivariate normal distribution manually

Question

I want to calculate the density of a multivariate normal distribution manually. As inputs of my function, I have x which is a n*p matrix of data points, a vector mu with n means and a covariance matrix sigma of dim p*p.

I wrote the following function for this:

`dmnorm <- function(mu, sigma, x){
k <- ncol(sigma)
x <- t(x)
dmn <- exp((-1/2)*t(x-mu)%*%solve(sigma)%*%(x- 
mu))/sqrt(((2*pi)^k)*det(sigma))  
return(dmn)
}`

My own function gives me a matrix of n*n. However, I should get a vector of length n.

In the end, I want the same results as I get from using the dmvnorm() function from the mvtnorm package. What's wrong with my code?

Answer 1

The expression t(x-mu)%*%solve(sigma)%*%(x- mu) is p x p, so that's why your result is that size. You want the diagonal of that matrix, which you can get using

diag(t(x-mu)%*%solve(sigma)%*%(x-mu))

so the full function should be

dmnorm <- function(mu, sigma, x){
  k <- ncol(sigma)
  x <- t(x)
  dmn <- exp((-1/2)*diag(t(x-mu)%*%solve(sigma)%*%(x- 
    mu)))/sqrt(((2*pi)^k)*det(sigma))  
  dmn
}