Integrate out (covariance) matrix in Normal-Wishart distribution

Question

In Gelman's Bayesian Data Analysis Chapter 3.6, he introduces the multivariate normal with unknown mean and variance, with the priors

$\Sigma\sim \text{Inv-Wishart}_{\nu_0}(\Lambda_0^{-1})$

$\mu\rvert \Sigma \sim N(\mu_0, \Sigma/\kappa_0)$

and the joint prior

$p(\mu, \Sigma)\propto \lvert \Sigma\rvert^{-((\nu_0+d)/2+1)}\exp\Big(-\frac{1}{2}\text{tr}(\Lambda_0\Sigma^{-1})-\frac{\kappa_0}{2}(\mu-\mu_0)'\Sigma^{-1}(\mu-\mu_0)\Big)$

Gelman then states that the (posterior) marginal is multivariate $t$ without elaboration.

In a paper that I found (A Bayes Approach for Combining Correlated Estimates - Seymour Geisser), Geisser had the following steps, also without elaboration (notations, parameterisation and context are different from Gelman, but a similar integration step is involved).

$p(\mu, \Sigma^{-1})\propto \Big\lvert \Sigma^{-1}\Big\rvert ^{(n-d-1)/2}\exp\Big(-\frac{1}{2}\text{tr}\Big(\Sigma^{-1}\sum_{j=1}^n (\mathbf{x}_j-\mu\mathbf{e})(\mathbf{x}_j-\mu\mathbf{e})'\Big)\Big)$

Integrating out $\Sigma^{-1}$ yields

$p(\mu)\propto \Big\lvert \sum_{j=1}^n(\mathbf{x}_j-\mu\mathbf{e})(\mathbf{x}_j-\mu\mathbf{e})'\Big\rvert ^{-n/2}$

The result can be further massaged to yield the $t$ density, but my question is specifically on the integration step only. The result $p(\mu)$ looks simple but I'm not sure how it's carried out.

Question: How do I obtain the marginal distribution of $p(\mu)$ by integrating out $\Sigma$? Any references to similar matrix integrations will be great. Thank you in advance!

Answer 1

First, notice that, if you use some properties of the trace operator,

\begin{align*} p(\mu, \Sigma) &\propto \lvert \Sigma\rvert^{-((\nu_0+d)/2+1)}\exp\Big(-\frac{1}{2}\text{tr}(\Lambda_0\Sigma^{-1})-\frac{\kappa_0}{2}(\mu-\mu_0)'\Sigma^{-1}(\mu-\mu_0)\Big) \\ &= \lvert \Sigma\rvert^{-((\nu_0+d)/2+1)}\exp\Big(-\frac{1}{2}\text{tr}(\Lambda_0\Sigma^{-1})-\text{tr}\left(\frac{\kappa_0}{2}(\mu-\mu_0)(\mu-\mu_0)'\Sigma^{-1}\right)\Big) \\ &= \lvert \Sigma\rvert^{-((\nu_0+d)/2+1)}\exp\Big(-\frac{1}{2}\text{tr}\left\{\Lambda_0 + \kappa_0(\mu-\mu_0)(\mu-\mu_0)' \right\} \Sigma^{-1}\Big) \\ &= \lvert \Sigma\rvert^{-\left(\frac{[\nu_0+1] + d+1}{2}\right)}\exp\Big[-\frac{1}{2}\text{tr}\left( \Psi \Sigma^{-1} \right)\Big] \\ &= \lvert \Sigma\rvert^{-\left(\frac{\nu_1 + d+1}{2}\right)}\exp\Big[-\frac{1}{2}\text{tr}\left( \Psi \Sigma^{-1} \right)\Big] \end{align*}

if you set $\Psi$ matrix to $\Lambda_0 + \kappa_0(\mu-\mu_0)(\mu-\mu_0)' $ and $\nu_1 = \nu_0 + 1$. So if you integrate with respect to $d \Sigma$, which is kind of misleading because all of the elements aren't unique, you'll get

$$ p(\mu) \propto \int \lvert \Sigma\rvert^{-\left(\frac{\nu_1 + d+1}{2}\right)}\exp\Big[-\frac{1}{2}\text{tr}\left( \Psi \Sigma^{-1} \right)\Big] d \Sigma = \frac{\Gamma_p(\nu_1/2) 2^{\nu_1p/2} }{\text{det}[\Psi]^{\nu_1/2}}. $$

Technically, $d \Sigma$ is short hand for integrating with respect to the diagonal and the lower-half elements of $\Sigma$. This is because you're integrating over the space of positive definite, symmetric matrices. I can skip a lot of the details because I am just recognizing the Inverse-Wishart density and then looking up the normalizing constant on Wikipedia.

Edit:

Just re-reading your question now and it looks like you're more interested in the marginal posterior. This is the marginal prior. You can use the same move, but I'll add the details in a little bit.

Edit Number 2:

There are a few extra difficulties that I wanted to mention with the recognition of the $t$ density. This site doesn't have many answered questions like this, so I figured I could do it all out.

\begin{align*} p(\mu) &\propto \text{det}\left(\Lambda_0 + \kappa_0(\mu-\mu_0)(\mu-\mu_0)'\right)^{-\frac{\nu_1}{2}} \\ &= \text{det}[\Lambda_0]^{-\frac{\nu_1}{2}} \text{det}\left(1 + \kappa_0(\mu-\mu_0)' \Lambda_0^{-1} (\mu-\mu_0) \right)^{-\frac{\nu_1}{2}} \tag{*} \\ &\propto \text{det}\left(1 + \kappa_0(\mu-\mu_0)' \Lambda_0^{-1} (\mu-\mu_0) \right)^{-\frac{\nu_1}{2}} \\ &= \text{det}\left(1 + \frac{1}{\nu_1 - p}(\mu-\mu_0)' \left[\frac{1}{[\nu_1 - p] \kappa_0}\Lambda_0\right]^{-1} (\mu-\mu_0) \right)^{-\frac{[\nu_1-p]+p}{2}} \\ &\propto t_{\nu_1 - p}\left(\mu_0, \frac{1}{[\nu_1 - p] \kappa_0}\Lambda_0\right). \end{align*} where the line marked by the * follows by the matrix determinant lemma.

You might think that it's weird you lost a few degrees of freedom, but you can convince yourself it's true if you check the variance is what it should be. The variance of a $t_{\nu_1 - p}\left(\mu_0, \frac{1}{[\nu_1 - p] \kappa_0}\Lambda_0\right)$ is $$ \frac{1}{[\nu_1 - p] \kappa_0}\Lambda_0 \frac{\nu_1 - p}{\nu_1 - p-2} = \frac{1}{[\nu_1 - p-2] \kappa_0}\Lambda_0 $$ which is the same as the one we find with the law of total variance: $$ E\left(\text{Var}\left[\mu \mid \Sigma \right]\right) = E[\Sigma]/\kappa_0 = \Lambda_0 \frac{1}{\kappa_0[\nu_0 - p - 1]}. $$