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.
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
.