I am trying to obtain the eigenvalue distribution of a finite dimensional Wishart matrix. Let $A_{n\times n}\sim\mathbb{W}(\Sigma_{n\times n},m)$ where $\mathbb{W}(\Sigma_{n\times n},m)$ denotes the Wishart distribution with covariance $\Sigma_{n\times n}$ and degrees of freedom $m$. I know that the result in the asymptotic case as $n,m\to\infty$ and $n/m\to c$ is given by the Marchenko–Pastur (M-P) density. My research led me to the following book:

T. W. Anderson,

"An Introduction to Multivariate Statistical Analysis", Ed. 2, John Wiley & Sons, Inc., 1984

It is my understanding that this is a classical textbook on multivariate statistics and a lot of its results are very well cited. One that was of particular interest to me is in page 534, Theorem 13.3.2 (notation changed to match the usage above):

If $A\ (n\times n)$ has the distribution $\mathbb{W}(I,m)$, then the characteristic roots $(l_1\geq l_2\geq\ldots\geq l_n\geq 0)$ have the [following] density over the range when the density is not 0.

$$\frac{\pi^{n^2/2}\prod_{i=1}^n l_i^{(m-n-1)/2}\exp\left(-\frac{1}{2}\sum_{i=1}^n l_i\right)\prod_{i\lt j}(l_i-l_j)}{2^{nm/2}\Gamma_n(m/2)\Gamma_n(n/2)}$$

where $\Gamma_n(x)$ is the multivariate Gamma function defined as

$$\Gamma_n(x)=\pi^{n(n-1)/4}\prod_{i=1}^n \Gamma\left(x+(1-i)/2\right)$$

I have the following questions:

- The expression for the density is a bit strange to me. Am I right in understanding that if there is an underlying continuous density $\rho(x)$, the expression above gives the density at a particular point $l_j$ as $\rho(l_j)$? If so, how do I flesh out $\rho(x)$ from the complicated expression?
- It seems to me that the expression requires knowledge of the eigenvalues $l_1,\ldots, l_n$ to calculate the density, and as such is less useful than a function $\rho(x)$ (like the M-P density), which is free of the individual eigenvalues. For example, if I needed to know the density for $n=10000$, I'll have to first evaluate the eigenvalues for a $10000\times 10000$ matrix before using the expression. In that case, bin counting or histogram would be simpler to use than this expression.
- Coming to evaluating the above density numerically, the denominator seems to be humungous ($10^{5000}$ish) when compared to the numerator ($10^{300}$ish) for nominal values of say, $n=30$ and $m=100$. As such, the density evaluates to
*almost*0 everywhere. I haven't posted the code as I believe it is my interpretation of the expression above (questions 1 & 2) that are incorrect, and not my implementation. However, I can provide code in Mathematica to evaluate the above expression, if anyone wants it.