2
$\begingroup$

My RandomVariate fails to evaluate when sampling from a multinormal distribution, any idea what's happening here?

This is Mathematica 12.1 on Mac:

d = 6;
sigma = {{0.2456589529728737`, 0.06467112508291382`, 
    0.07381608109139866`, 0.01761648518373636`, 0.05925859721052876`, 
    0.045898960400932734`}, {0.06467112508291387`, 
    0.35501955566599797`, 
    0.05242856376555134`, -0.009970312195220573`, 
    0.015359942596352098`, 
    0.02448314785173377`}, {0.07381608109139871`, 
    0.05242856376555134`, 0.3026715722222353`, -0.0819524814567561`, 
    0.07919970253622352`, -0.033627335100296744`}, \
{0.01761648518373638`, -0.009970312195220557`, -0.08195248145675606`, 
    0.5519153257149109`, -0.1068423071796184`, 
    0.35028088018797104`}, {0.05925859721052879`, 
    0.015359942596352072`, 
    0.07919970253622348`, -0.10684230717961839`, 0.29853333010040306`,
     0.002543954098042138`}, {0.0458989604009327`, 
    0.02448314785173373`, -0.033627335100296814`, 
    0.35028088018797104`, 0.002543954098042013`, 0.6962012633235793`}};
ii = IdentityMatrix[d];
RandomVariate@MultinormalDistribution[sigma] (* fails *)
RandomVariate@MultinormalDistribution[sigma + ii] (* works *)

The example above was generated by code below, it's a random covariance matrix with geometrically decaying eigenvalues, so not even that close to singular --

randomRotation[n_] := Module[{M, z, q, r, d, ph},
   z = RandomVariate[NormalDistribution[0, 1], {n, n}];
   {q, r} = QRDecomposition[z];
   d = Diagonal[r];
   ph = d/Abs[d];
   M = q*ph; (* Note, determinant may be -1 giving "pseudo-rotation", 
   switch 2 rows of the matrix to guarantee true rotation *)
   indices = If[Det[M] > 0, Range[n], {2, 1}~Join~Range[3, n]];
   M[[indices]]
   ];

d = 20;
rot = randomRotation[d];
ii = IdentityMatrix[d];

sigma = rot.DiagonalMatrix[Table[1/i, {i, 1, d}]].Inverse[rot];
RandomVariate@MultinormalDistribution[sigma] (* fails *)
RandomVariate@MultinormalDistribution[sigma+ii] (* works *)
```
$\endgroup$
  • 2
    $\begingroup$ "...symmetric positive definite..." for the covariance. Sigma does not meet the condition. That said, an error should be for that, instead it seems to fail in the definition of the distribution itself (that is, MultinormalDistribution[sigma] should note the condition failure). You could try something like RandomVariate@MultinormalDistribution[N[Rationalize[sigma, 0], 8]] $\endgroup$ – ciao Oct 27 at 22:43
  • $\begingroup$ ah, thanks...it looks like symmetry failure from numerical noise. On Mac the error was not very informative, it just said MultinormalDistribution called with 1 argument; 2 arguments expected $\endgroup$ – Yaroslav Bulatov Oct 27 at 22:58
  • $\begingroup$ For whatever it's worth on Windows one gets the same 1 vs 2 arguments error. Also, see the answer to mathematica.stackexchange.com/questions/152987/… which gives a function to force the approximately symmetric matrix to be symmetric. $\endgroup$ – JimB Oct 28 at 0:27

Your Answer

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

Browse other questions tagged or ask your own question.