# Semi-conjugate inverse Wishart posterior, can we obtain the marginal?

In Hoff's text (A First Course in Bayesian Statistical Methods), he uses a semi-conjugate inverse-Wishart prior for the covariance matrix of a multivariate normal process. In equation 7.9, he has the following posterior for the covariance matrix $$\Sigma$$:

$$\{\Sigma\vert \mathbf{y}_1, ..., \mathbf{y}_n, \mathbf{\theta}\}\sim \text{inverse-Wishart}\Big(\nu_0+n, [\mathbf{S}_0+\mathbf{S}_\theta]^{-1}\big)$$,

where $$\mathbf{S}_0$$ is the prior scale matrix and $$\mathbf{S}_\theta=\sum_{i=1}^n (\mathbf{y}_i-\theta)(\mathbf{y}_i-\theta)'$$. If I'm not wrong, this posterior can also be viewed as the posterior for $$\Sigma$$ given a diffuse prior for $$\theta$$ and an inverse Wishart prior for $$\Sigma$$.

My question is this: Since the posterior here depends on $$\theta$$, is it possible to remove this conditioning by integrating out $$\theta$$ from the posterior inverse-Wishart distribution? This would yield some expression for $$\{\Sigma\vert \mathbf{y}_1, ..., \mathbf{y}_n\}$$, or the marginal posterior for $$\Sigma$$ if we assume a diffuse prior for $$\theta$$.

I don't think this yields a closed form solution, is this the case?

No, you cannot. Since in your problem setting, $$\Sigma | \vec{\mathbf{y}}$$ and $$\theta | \vec{\mathbf{y}}$$ are by no means independent. Thus we cannot simply marginalize it.

Here the problem setting is follows:

\begin{aligned} \mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma &\sim \mathcal{N}\left( \theta, \Sigma \right) \\[10pt] \Sigma &\sim \mathcal{W}^{-1} \left( \nu_0, \mathbf{S}_0 \right) \end{aligned}

where $$\mathcal{W}^{-1}(\nu, \mathcal{\Psi})$$ denotes the inverse-Wishart distribution, which is a generalization of inverse gamma$$(\alpha, \beta)$$ distribution when dimension $$p=1$$, with $$\nu = 2\alpha, \mathcal{\Psi} = 2\beta$$.

Note that we have assumed a given (nonrandom) $$\theta$$. This is why we can write

$$\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma \sim \mathcal{N}\left( \theta, \Sigma \right)$$

in the first place, rather than $$\mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma, \theta \sim \mathcal{N}\left( \theta, \Sigma \right)$$.

## Posterior when both $$\theta, \Sigma$$ : unknown

Alternatively, let us suppose $$\theta$$ is random, so that we can assume a prior on $$\theta$$. This is exactly where your idea may kick in. Note that when we start assuming

\begin{aligned} \mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma, \theta &\sim \mathcal{N}\left( \theta, \Sigma \right) \\[10pt] \Sigma &\sim \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \end{aligned}

we also simply get the posterior of $$\Sigma | \vec{\mathbf{y}}, \theta$$ , since by the Bayes rule,

\begin{aligned} p(\Sigma | \vec{\mathbf{y}}, \theta) &= \frac{p(\vec{\mathbf{y}}|\Sigma, \theta) p(\Sigma | \theta)}{\int_{\tilde{\Sigma}} p(\vec{\mathbf{y}}|\tilde{\Sigma}, \theta) p(\tilde{\Sigma}| \theta)} \\[8pt] &\propto \mathcal{N} \left( \nu_0 + n, [\mathbf{S} + \mathbf{S}_\theta]^{-1} \right) \end{aligned}

as in the equation (7.9). Note that unlike the above, we cannot say $$\Sigma|\theta \equiv \Sigma$$. This is exactly why your idea might need caveats. Rather, if we allow an additional information to be

$$\theta | \Sigma \sim \mathcal{N} (\theta_0, \Sigma/m)$$

for a nonrandom hyperparameters $$\theta_0$$,$$m$$, things get simpler, i.e. we can get a closed-form solution.

## Posterior when $$\theta | \Sigma$$ : normal

The joint density of $$\theta$$, $$\Sigma$$ is commonly referred to as Normal-inverse-Wishart distribution, which is simply formulated as

$$\theta, \Sigma \sim \mathcal{N} (\theta_0, \Sigma/m) \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \equiv NIW(\theta, m, \mathbf{S}_0, \nu_0)$$

It can be shown that the posterior distribution is also the Normal-inverse-Wishart, as

\begin{aligned} \theta | \Sigma , \vec{\mathbf{y}} &\sim \mathcal{N} (\theta_n, \Sigma/m_n) \\[8pt] \Sigma | \vec{\mathbf{y}} &\sim \mathcal{W}^{-1} (\nu_n, \mathbf{S}_n) \end{aligned}

where

\begin{aligned} \theta_n &= \frac{m \theta_0 n + \bar{\mathbf{y}} }{ m_n } \\ m_n &= m + n\\ \nu_n &= \nu_0 + n \\ \mathbf{S}_n &= \mathbf{S}_0 + \mathcal{S} + \frac{mn}{m+n} (\bar{\mathbf{y}} - \theta_0) (\bar{\mathbf{y}} - \theta_0)^T \end{aligned}

where $$\mathcal{S} := (\vec{\mathbf{y}}-\bar{\mathbf{y}}\mathbb{1})(\vec{\mathbf{y}}-\bar{\mathbf{y}}\mathbb{1})^T$$.

So you are now assuming the model

\begin{aligned} \mathbf{y}_1, \cdots, \mathbf{y}_n | \Sigma, \theta &\sim \mathcal{N}\left( \theta, \Sigma \right) \\[10pt] \Sigma &\sim \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \\[10pt] \theta &\propto 1 \\[10pt] \end{aligned}

where definitely the improper prior $$\theta$$ on $$p$$-dimensional real space $$\mathbb{R}^p$$ assumed. Most Bayesians accept improper priors if the resulting posterior is proper(see here).

So let us see if this setting leads to the proper posterior. We have

\begin{aligned} p(\theta, \Sigma) &\propto \mathcal{W}^{-1} (\nu_0, \mathbf{S}_0) \prod_{i=1}^n \mathcal{N}_\mathbf{y}\left( \theta, \Sigma \right) \\[10pt] \end{aligned}

thus we have

$$\Sigma , \theta | \vec{\mathbf{y}} \sim \mathcal{W}^{-1} (\nu_0 + n, [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1})$$

Note that the difference between here and in (7.9) is the LHS in (7.9) is $$\Sigma | \theta, \vec{\mathbf{y}}$$, whereas we have joint here. Thus, we can integrate out $$\theta$$, (kind of) as what you said,

\begin{aligned} p(\Sigma | \vec{\mathbf{y}}) &= \int_\Theta \mathcal{W}^{-1} (\nu_0 + n, [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1}) d\theta \\[10pt] &= \int_{\mathbb{R}^p} \frac{\left|\mathbf{S}_0 + \mathbf{S}_\theta \right|^{-(\nu_0 + n)/2}}{2^{(\nu_0 + n) p/2}\Gamma_p(\frac{\nu_0 + n}{2})} \left| \Sigma \right|^{-(\nu_0 + n+p+1)/2} e^{-\frac{1}{2}\operatorname{tr}( [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1} \Sigma^{-1})} d\theta \\[10pt] &= \frac{\left| \Sigma \right|^{-(\nu_0 + n+p+1)/2}}{2^{(\nu_0 + n) p/2}\Gamma_p(\frac{\nu_0 + n}{2})} \int_{\mathbb{R}^p} \left|\mathbf{S}_0 + \mathbf{S}_\theta \right|^{-(\nu_0 + n)/2} e^{-\frac{1}{2}\operatorname{tr}( [\mathbf{S}_0 + \mathbf{S}_\theta]^{-1} \Sigma^{-1})} d\theta \end{aligned}

This integral does not converge, which means it is an improper posterior. Thus, most Bayesians may not accept this as the posterior.

To see why the integral does not converge in a little bit heuristic way, say for the one dimensional case($$p=1$$), the integral is the form

$$p(\sigma^2 | \vec{\mathbf{y}}) = C \cdot \int_\mathbf{R} \theta^{-(\nu_0 + n)} exp(-\theta^2) d\theta$$

for some constant $$C$$. It is known that the above integral converges iff

$$\nu_0 + n < 1$$

which is not attainable for $$n \ge 1$$ cases.

PS. Interestingly, we can show that the marginal prior density of the $$\theta$$ is multivariate $$t$$ distribution (reference: Bernardo and Smith[1]), where

$$\theta \sim T(\theta_0, \kappa, q)$$

where $$\kappa=\nu_0 \mathbf{S}_0/m, q=2\nu_0 - p + 1$$.

### Reference

[1] Bernardo and Smith, Bayesian theory(1994), p435

• Thanks for the detailed response! I have some clarifications regarding your section "Posterior when both are unknown". If we allow the additional information on theta to be improper diffuse so that p(theta)=1 instead of the standard multivariate normal, we yield the same posterior as equation 7.9. Conceptually, this is acceptable right? In this case, it would also be conceptually acceptable to think about integrating out theta to yield a marginal, and the only problem is that there is no closed form solutions? Aug 23 '19 at 1:19
• @bayes No, the real problem is, that does not yield a proper posterior. Please check the latest edit. Aug 23 '19 at 5:39