Solved – Multivariate normal log-likelihood computation in R

I'm trying to use the dmvnorm function from the mvtnorm package to compute the log-likelihood of a random normal vector. However, when I check the output of dmvnorm against my computation of the log-likelihood function the dmvnorm output is always almost exactly my result multiplied by 1.7. What's going on? Am I computing the likelihood wrong? Am I using the dmvnrom function wrong? Here's my code

x <- rmvnorm(1, sigma = C1) loglik <- dmvnorm(x, sigma = C1, log = T)  y <- t(x) deter2pi <- function(mat){     determinant(2*pi*mat)[[1]][1] } ll <- (-0.5*t(y)%*%solve(C1,y)-0.5*log(deter2pi(C1)))[1] loglik/ll # = 1.7 

And my covariance matrix

> round(C1,3)        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]  [1,] 1.000 0.743 0.545 0.421 0.309 0.215 0.159 0.120 0.084 0.059  [2,] 0.743 1.000 0.733 0.564 0.413 0.289 0.213 0.161 0.113 0.079  [3,] 0.545 0.733 1.000 0.766 0.561 0.393 0.290 0.220 0.154 0.108  [4,] 0.421 0.564 0.766 1.000 0.732 0.511 0.377 0.286 0.200 0.141  [5,] 0.309 0.413 0.561 0.732 1.000 0.695 0.515 0.390 0.273 0.192  [6,] 0.215 0.289 0.393 0.511 0.695 1.000 0.736 0.559 0.391 0.275  [7,] 0.159 0.213 0.290 0.377 0.515 0.736 1.000 0.758 0.531 0.372  [8,] 0.120 0.161 0.220 0.286 0.390 0.559 0.758 1.000 0.700 0.491  [9,] 0.084 0.113 0.154 0.200 0.273 0.391 0.531 0.700 1.000 0.699 [10,] 0.059 0.079 0.108 0.141 0.192 0.275 0.372 0.491 0.699 1.000 

By default, determinant returns the log of the determinant.

   x <- rmvnorm(1, sigma = C1)     loglik <- dmvnorm(x, sigma = C1, log = T)      y <- t(x)     deter2pi <- function(mat){         determinant(2*pi*mat,log=FALSE)[[1]][1]     }     ll <- (-0.5*t(y)%*%solve(C1,y)-0.5*log(deter2pi(C1)))[1]     loglik/ll 

does return the value 1 for the ratio.

Similar Posts:

Rate this post

Leave a Comment