If $S\sim \mathcal{W}_p(\Sigma, n)$, what is the distribution of the diagonal part of $S$, $\text{diag}(S)$?
Principal submatrices of Wishart distributed matrices are also Wishart distributed, implying that each diagonal element is marginally gamma distributed. But what is the joint distribution of these diagonal terms?
A refresher on the Wishart distribution:
Let $X_i \sim \mathcal{N}_p(0, \Sigma)$ be a vector independently drawn from p-variate normal distribution with mean $0$ and covariance $\Sigma$ for $i\in\{1,\ldots,p\}$. Then
$$ S \equiv \sum_{i=1}^{n}X_iX_i^T \sim \mathcal{W}_p(\Sigma, n)\,. $$
I have been able to determine the joint moment generating function (MFG) of diag($\Sigma$), and I will include the derivation here. I am a bit stuck at this point however, so feel free to skip to the bottom or ignore this work entirely if you think there is a better approach.
The meat of my work so far:
Let's denote the diagonal elements of $S$ as $\sigma \equiv \text{diag}(\Sigma)$, with $\sigma_i = S_{ii}$. Then we have that
$$ \sigma_i = \left(\sum_{j=1}^{n}X_jX_j^T\right)_{ii} = \sum_{j=1}^{n}(X_j)_i^2\,. $$
In vector notation, we have $$ \sigma = \sum_{i=1}^{n} X_i\odot X_i\,, $$ where $\odot$is the Hadamard (elementwise) product.
Let's consider the MGF of $\sigma$: $$ \phi_\sigma^n(t) \equiv \left\langle e^{t^T\sigma} \right\rangle_{p(\sigma)} = \left\langle e^{t^T\left(\sum_{i=1}^nX_i\odot X_i\right)} \right\rangle_{p(X_1, \ldots, X_n)} = \left\langle e^{t^T\left(X\odot X\right)} \right\rangle_{p(X)}^n\,, $$
where the first equality is given by the law of the unconscious statistician and the second is given by the independence of $X_i$ and $X_j$ for $i \neq j$.
Considering just $$ \phi(t) \equiv \phi_\sigma^1(t) = \left\langle e^{t^T\left(X\odot X\right)} \right\rangle_{p(X)}\,, $$ we have $$ \begin{align} \phi(t) &= \int_{-\infty}^{\infty} dX \ p(X) e^{t^T\left(X\odot X\right)} \\ &= \int_{-\infty}^{\infty} dX \ \frac{1}{\sqrt{2\pi}} |\Sigma|^{-1/2} \exp\left(-\frac{1}{2} X^T\Sigma^{-1}X\right) \exp\left(t^T\left(X\odot X\right)\right)\,. \end{align} $$
If we let $D \equiv \begin{bmatrix}t_1 & &0\\ &\ddots& \\ 0&&t_p \end{bmatrix}$, then we can write $t^T\left(X\odot X\right) = X^T D X$:
$$ \begin{align} \phi(t) &= \int_{-\infty}^{\infty} dX \ \frac{1}{\sqrt{2\pi}} |\Sigma|^{-1/2} \exp\left(-\frac{1}{2} X^T\Sigma^{-1}X\right) \exp\left(X^T D X\right) \\ &= \int_{-\infty}^{\infty} dX \ \frac{1}{\sqrt{2\pi}} |\Sigma|^{-1/2} \exp\left(-\frac{1}{2} X^T(\Sigma^{-1} - 2D)X\right)\,. \end{align} $$
If $\Sigma$ is positive definite, we have a finite radius $\epsilon > 0$ such that $\Sigma^{-1} - 2D$ is positive definite for $|t_i| < \epsilon \ \forall i$, and we can find a full-rank $p\times p$ matrix $B$ such that $B^TB = \Sigma^{-1} - 2D$ in this radius.
Let's change variables, letting $Y \equiv BX$. Including the Jacobian, we have $$ \begin{align} \phi(t) &= \int_{-\infty}^{\infty} dX \ \frac{1}{\sqrt{2\pi}} |\Sigma|^{-1/2} \exp\left(-\frac{1}{2} X^T(\Sigma^{-1} - 2D)X\right) \\ &= |\Sigma|^{-1/2} |B|^{-1} \int_{-\infty}^{\infty} dY \ \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{1}{2} Y^TY\right) \\ &= |\Sigma|^{-1/2} |B|^{-1} \\ &= |\Sigma|^{-1/2} |\Sigma^{-1} - 2D|^{-1/2} \\ &= |I_p - 2\Sigma D|^{-1/2}\,. \end{align} $$
Finally, we have (TLDR): $$ \phi_\sigma^n(t) = |I_p - 2\Sigma D|^{-\frac{n}{2}}\,. $$
This feels like a neat little result in of itself. I am fairly confident that this is correct, since it reduces to the MGF of the gamma distribution for $p=1$, and it nicely reduces to the product of several gamma MFGs if $\Sigma$ is diagonal.
However, I don't really know where to go from here. I can't really compute the inverse Laplace transform by inspection here, and I quickly get very lost when trying to do it explicitly.
Is there a known result anybody can point me to? Is there a nice method of retrieving a joint density for $\sigma$ from this expression? Thank you very much for the help!