1
$\begingroup$

I'm drawing from the Inverse-Wishart Distribution, but I got the following error message:

InverseWishartMatrixDistribution::posdefprm: 
   The value {{768.252, 317.247, 737.618}, {317.247, 180.387, 305.871}, 
     {737.618, 305.871, 708.24}} at position 2 in 
    InverseWishartMatrixDistribution[10, 
     {{768.252, 317.247, 737.618}, {317.247, 180.387, 305.871}, 
      {737.618, 305.871, 708.24}}] is expected to be a symmetric positive
     definite matrix.

Well, when I do SymmetricMatrixQ and PositiveDefiniteMatrixQ, I get True for the matrix shown... It's probably that Mathematica has approximated the true values in the error message.

The problem lies, I suspect elsewhere. The code I use is the following:

RandomVariate[
 InverseWishartMatrixDistribution[nu, mean*(nu + Length[mean] + 1)]]

If mean is PD and Sym, then mean*(nu + Length[mean] + 1) with nu=10 and Length[mean]=3 should also be...

I make sure that mean is PD and Sym with the another programme, let us call it SPDMatrix.

So, I think that the problem maybe that InverseWishartMatrixDistribution may somehow use greater machine precision functions like SymmetricMatrixQ and PositiveDefiniteQ than those available to the users...

I've managed to replicate the problem in a simpler way:

In[102]:= RandomSeed[1234];
list = RandomReal[{1, 3}, {100000, 3, 3}];

In[105]:= 
list2 = ParallelTable[
    SPDMatrix[list[[i]]], {i, 1, Length[list]}]*14;

In[107]:= Map[SymmetricMatrixQ, list2] // Total

Out[107]= 100000 True

In[106]:= Map[PositiveDefiniteMatrixQ, list2] // Total

Out[106]= 22234 False + 77766 True

In[108]:= 
list3 = ParallelTable[
   SPDMatrix[list[[i]]], {i, 1, Length[list]}];

In[109]:= Map[PositiveDefiniteMatrixQ, list3] // Total

Out[109]= 100000 True

The difference between list2 and list3 is that in list 2 I've multiplied by 14... This must interfere with the machine precision.

$\endgroup$
  • $\begingroup$ Can you give an explicit SPD matrix mat such that RandomVariate[InverseWishartMatrixDistribution[n, mat]] throws an error message? Copying your matrix {{768.252, 317.247, 737.618}, {317.247, 180.387, 305.871}, {737.618, 305.871, 708.24}} into Mathematica and using it in the distribution shows no problems. $\endgroup$ – J. M. will be back soon Apr 1 at 9:20
  • $\begingroup$ @J.M.isslightlypensive I'm running my code on a server, to which I have no interface. (I'm windows, and the server is linux, so I just run a .m file). So, the only way I have of detecting this is by using SymmetricMatrixQ, and PositiveDefiniteQ, which do NOT detect the problem.... $\endgroup$ – An old man in the sea. Apr 1 at 9:23
  • $\begingroup$ Can't you modify your code such that it returns InputForm[mean (nu + Length[mean] + 1)] whenever a failure is encountered? $\endgroup$ – J. M. will be back soon Apr 1 at 9:24
  • $\begingroup$ @J.M.isslightlypensive Probably. Do you know of a way to create a condition for when Mathematica emits an error message to print that Inputform? Otherwise, it would need to print possibly a few thousands of matrices before we get the one we want... $\endgroup$ – An old man in the sea. Apr 1 at 9:32
  • $\begingroup$ Sure, insert a line like Check[RandomVariate[InverseWishartMatrixDistribution[nu, mean*(nu + Length[mean] + 1)]], Return[InputForm[mean*(nu + Length[mean] + 1)]], General::posdefprm] in the relevant section of your code. (Thankfully, you were thoughtful enough to mention what error message was being returned.) $\endgroup$ – J. M. will be back soon Apr 1 at 10:43
2
$\begingroup$

So, instead of just making sure that mean is SPD, I could just create a new function to draw from Inverse-Wishart, but making sure that mean*(nu + Length[mean] + 1) is SPD... Here's an example

drawIW[mean_, s2_] := Module[{matrix},
   nu = 10;
   matrix = SPDMatrix[mean*(nu + Length[mean] + 1)];
   (*We need this, since sometimes multiplying by nu+Length[
   mean]+1 will return a non-
   SPD matrix due to machine precision considerations... *)

   RandomVariate[InverseWishartMatrixDistribution[nu, matrix]]
   ];
$\endgroup$

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.