} \Gamma_p \left(\frac{\nu}{2} \right)}
           \exp\left( -tr(\Sigma S^{-1}) / 2 \right)

If :math:`S \sim W_p^{-1}(\nu, \Psi)` (inverse Wishart) then
:math:`S^{-1} \sim W_p(\nu, \Psi^{-1})` (Wishart).

If the scale matrix is 1-dimensional and equal to one, then the inverse
Wishart distribution :math:`W_1(\nu, 1)` collapses to the
inverse Gamma distribution with parameters shape = :math:`\frac{\nu}{2}`
and scale = :math:`\frac{1}{2}`.

Instead of inverting a randomly generated Wishart matrix as described in [2],
here the algorithm in [4] is used to directly generate a random inverse-Wishart
matrix without inversion.

.. versionadded:: 0.16.0

References
----------
.. [1] M.L. Eaton, "Multivariate Statistics: A Vector Space Approach",
       Wiley, 1983.
.. [2] M.C. Jones, "Generating Inverse Wishart Matrices", Communications
       in Statistics - Simulation and Computation, vol. 14.2, pp.511-514,
       1985.
.. [3] Gupta, M. and Srivastava, S. "Parametric Bayesian Estimation of
       Differential Entropy and Relative Entropy". Entropy 12, 818 - 843.
       2010.
.. [4] S.D. Axen, "Efficiently generating inverse-Wishart matrices and
       their Cholesky factors", :arXiv:`2310.15884v1`. 2023.

Examples
--------
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.stats import invwishart, invgamma
>>> x = np.linspace(0.01, 1, 100)
>>> iw = invwishart.pdf(x, df=6, scale=1)
>>> iw[:3]
array([  1.20546865e-15,   5.42497807e-06,   4.45813929e-03])
>>> ig = invgamma.pdf(x, 6/2., scale=1./2)
>>> ig[:3]
array([  1.20546865e-15,   5.42497807e-06,   4.45813929e-03])
>>> plt.plot(x, iw)
>>> plt.show()

The input quantiles can be any shape of array, as long as the last
axis labels the components.

Alternatively, the object may be called (as a function) to fix the degrees
of freedom and scale parameters, returning a "frozen" inverse Wishart
random variable:

>>> rv = invwishart(df=1, scale=1)
>>> # Frozen object with the same methods but holding the given
>>> # degrees of freedom and scale fixed.

c