Computing the MultinormalDistribution density gives sometimes different values from in-built function

Question

I'm trying to define the density of the Multivariate Normal distribution, since it will be faster to compute for greater dimensions (vectors of size 50 or more), specially when I have to do those computations many times (10,000).

I've defined the density as :

multinormalDens[x_, mean_, var_] := Module[{},
   Det[2*Pi*var]^(-0.5)*Exp[-0.5*(x - mean).Inverse[var].(x - mean)]
   ];

I use the Wishart to provide with several positive definite covariance matrices.

varmat = RandomVariate[
   WishartMatrixDistribution[52, IdentityMatrix[50]], 10000];
meanmat = Table[RandomReal[{0, 5}, 50], {i, 1, 10000}];

They will all be evaluated for the same point:

x0 = RandomReal[{0, 5}, 50];

Now running the comparison for 10,000 tries, I get:

In[43]:= Table[
  multinormalDens[x0, meanmat[[i]], varmat[[i]]] == 
   PDF[MultinormalDistribution[meanmat[[i]], varmat[[i]]], x0], {i, 1,
    Length[varmat]}] // Total

Out[43]= 6643 False + 3357 True

How come I have 6643 different outcomes? Are they from the matrix inversion?

Edit:

data = Table[
   multinormalDens[x0, meanmat[[i]], varmat[[i]]] - 
    PDF[MultinormalDistribution[meanmat[[i]], varmat[[i]]], x0], {i, 
    1, Length[varmat]}];



ListPlot[data]
ListPlot[data, PlotRange -> All]

Answer 1

This is less about accuracy (it is already good) but about getting rid of the General::munfl messages. Before the final number is computed, we should check whether the logarithm of the result is below (a multiple of) Log[$MinMachineNumber]; if yes, we simply set the result to 0.

multinormalDens[x_, mean_, var_] := Module[{t, det},
   t = 0.5*(x - mean).LinearSolve[var, (x - mean)];
   det = Det[2*Pi*var]^(0.5);
   If[-t - Log[det] < 10. Log[$MinMachineNumber], 0., Exp[-t]/det]
   ];