# Computing the Multivariate Normal 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. – b.gatessucks Feb 12 at 8:43
• @b.gatessucks I've plotted a graphic. ;) – An old man in the sea. Feb 12 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. – Henrik Schumacher Feb 12 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. – Henrik Schumacher Feb 12 at 8:53
• Your definition of meanmat can be simplified to meanmat = RandomReal[{0, 5}, {10000, 50}]; – Bob Hanlon Feb 12 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]