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

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]


• Maybe it is interesting to look at the distribution of the differences. Commented Feb 12, 2019 at 8:43
• @b.gatessucks I've plotted a graphic. ;) Commented Feb 12, 2019 at 8:49
• These (absolute) errors are really tiny. (The relative errors might be much larger...) In terms of probability, they may certainly be neglected. Commented Feb 12, 2019 at 8:51
• There are two major source for numerical errors: The matrix inversion is one. But more important are probably the underflows caused by Exp[-(...)] with large arguments. Commented Feb 12, 2019 at 8:53
• Your definition of meanmat can be simplified to meanmat = RandomReal[{0, 5}, {10000, 50}]; Commented Feb 12, 2019 at 15:19

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]