Drawing from InverseWishartMatrixDistribution. A machine precision error?

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.

• 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. – J. M. will be back soon Apr 1 at 9:20
• @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.... – An old man in the sea. Apr 1 at 9:23
• Can't you modify your code such that it returns InputForm[mean (nu + Length[mean] + 1)] whenever a failure is encountered? – J. M. will be back soon Apr 1 at 9:24
• @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... – An old man in the sea. Apr 1 at 9:32
• 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.) – J. M. will be back soon Apr 1 at 10:43

drawIW[mean_, s2_] := Module[{matrix},