Search code examples
rprobabilityevaluationnormal-distributionlog-likelihood

Efficiently evaluate Multivariate Normal


I want to evaluate datapoints that arise from multivariate normal densities. I have to evaluate each datapoint with respect to different means and covariance matrices. I have two means for each observation to evaluate the likelihood. Also, I have two different variance covariance matrices overall. For now, I am considering two-dimensional normal distributions only.

Basically, I have to do a lot of multivariate likelihood evaluations and I am looking for a way to do this faster. Here's some example code (data below):

N <- 10 #number of observations
G <- 2  #number of means per observation / variance - covariance matrices

ll <- array(NA, c(N,G)) #shell for the log likelihoods

for(ii in 1:N){ #loop over data-points
  for(gg in 1:G){ #loop over groups of means / var-cov matrices

    ll[ii,gg] <- mixtools::logdmvnorm(data[ii,], #evaluate data by observation
                               mu = means[[gg]][ii,],  #mean by group by obs.
                               sigma = Sigma[[gg]])    #var-cov matrix by group

  }
}

What I want to do is the following: Take the first data-point, evaluate it using mean A from observation 1 and covariance matrix A. Evaluate it with mean B from observation 1 and covariance matrix B. Take the second data-point, evaluate it with respect to mean A belonging to observation 2 / covariance matrix A. Then evaluate it with mean B from observation 2 / covariance matrix B and so on.

I have prepared 10 datapoints as well as 10*2 mean vectors and 2 variance covariance matrices here. It is not necessary to keep the list structure, it just arose naturally in the coding process.

In an univariate setting, it is possible to have sufficiently fast performance by using the fact that dnorm() is vectorized. Hence, N iterations are not necessary in this case.

Thanks!


Solution

  • Found a solution:

    mnormt::dmnorm is vectorized.