2
$\begingroup$

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]

enter image description here enter image description here

$\endgroup$
  • $\begingroup$ Maybe it is interesting to look at the distribution of the differences. $\endgroup$ – b.gates.you.know.what Feb 12 at 8:43
  • $\begingroup$ @b.gatessucks I've plotted a graphic. ;) $\endgroup$ – An old man in the sea. Feb 12 at 8:49
  • 1
    $\begingroup$ These (absolute) errors are really tiny. (The relative errors might be much larger...) In terms of probability, they may certainly be neglected. $\endgroup$ – Henrik Schumacher Feb 12 at 8:51
  • 1
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Feb 12 at 8:53
  • $\begingroup$ Your definition of meanmat can be simplified to meanmat = RandomReal[{0, 5}, {10000, 50}]; $\endgroup$ – Bob Hanlon Feb 12 at 15:19
3
$\begingroup$

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]
   ];
$\endgroup$
  • $\begingroup$ Henrik, thanks for the answer. What are the General::munfl messages ? They don't show up in my notebook... $\endgroup$ – An old man in the sea. Feb 12 at 9:41
  • 1
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Feb 12 at 10:13
  • $\begingroup$ Thanks Henrik ;) $\endgroup$ – An old man in the sea. Feb 12 at 10:20
  • $\begingroup$ Always at your service! =) $\endgroup$ – Henrik Schumacher Feb 12 at 10:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.