# 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.gates.you.know.what 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

## 1 Answer

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]
];

• Henrik, thanks for the answer. What are the General::munfl messages ? They don't show up in my notebook... – An old man in the sea. Feb 12 at 9:41
• Yes, it's likely that they do not appear if you use version 11.2 or older. Wolfram Research added these messages due to introduction of a different underflow handling in version 11.3. The aim is to be closer to the IEEE 754 standard for floating point arithmetic. – Henrik Schumacher Feb 12 at 10:13
• Thanks Henrik ;) – An old man in the sea. Feb 12 at 10:20
• Always at your service! =) – Henrik Schumacher Feb 12 at 10:21