Search code examples
wolfram-mathematicamathematica-8

difference in mathematica between Inverse[_] and (_)^(-1) for WishartDistribution


Does anyone know why the following random distributions of matrices generate different plots? (This is code to generate a plot of the PDFs for first cells from a set of 10x10 matrices sampled using an inverse Wishart distribution; amazingly, the plots are different depending on the way one performs the matrix inverse - and it seems the right plots are obtained by Inverse[_], why?)

base code:

   << MultivariateStatistics`;
   Module[{dist, p, k, data, samples, scale, graphics, distribution},
     p = 10;
     k = 13;
     samples = 500;
     dist = WishartDistribution[IdentityMatrix[p], k];
     (* a samples x p x p array *)
     data = Inverse[#] & /@ RandomVariate[dist, samples];
     (* distribution graphics *)
     distribution[i_, j_] := Module[{fiber, f, mean, rangeAll, colorHue},
       fiber = data[[All, i, j]];
       dist = SmoothKernelDistribution[fiber];
       f = PDF[dist];
       Plot[f[z], {z, -2, 2}, 
        PlotLabel -> ("Mean=" <> ToString[Mean[fiber]]), 
        PlotRange -> All]
       ];
     Grid @ Table[distribution[i, j], {i, 1, 3}, {j, 1, 5}]
     ]

code variant: above, change line

data = Inverse[#] & /@ RandomVariate[dist, samples];

by this

data = #^(-1) & /@ RandomVariate[dist, samples];

and you will see the plotted distributions are different.


Solution

  • Inverse computes a matrix inverse, i.e. if a is a square matrix, then Inverse[a].a will be the identity matrix.

    a^(-1) is the same as 1/a, i.e. it gives you the reciprocal of each matrix element. The ^ operator gives powers element-wise. If you want a matrix power, use MatrixPower.