Unsolvable Integral?

Is the following integral solvable?

$$P(X) = \int^{\infty}_{-\infty} \int^{\infty}_{-\infty} P(X|\mu,K)P(\mu|K)P(K) d\mu dK$$

with

$$P(K) = \frac{|K| ^{(v-d-1)/2}}{2^{vd/2}|V|^{v/2}\Gamma_d|\frac{v}{2}|} e^{-tr(V^{-1}K)/2}$$ (Wishart distribution)

$$P(\mu|K) = \frac{|\lambda_0K|^{1/2}}{(2\pi)^{d/2}}e^{-0.5([\mu - \mu_0]^T \lambda_0K[\mu-\mu_0])}$$

(Gaussian distribution)

$$P(X|\mu,K) = \frac{|K|^{1/2}}{(2\pi)^{d/2}}e^{-0.5([X- \mu]^T K[X-\mu])}$$

(Gaussian distribution - except not really - because $K$ and $\mu$ are random variables)

with $K$ being a matrix variable and $X$ and $\mu$ being vector variables

• This looks like the posterior predictive distribution for the conjugate Bayesian model of the multivariate Gaussian distribution. Is it that? – Stéphane Laurent Jun 30 '17 at 14:43
• in $P(K)$ is the $n$ supposed to be a $\nu$? – Taylor Jul 1 '17 at 2:33
• And you should only pick either $d$ or $p$ but not both. – Taylor Jul 1 '17 at 2:43
• @Taylor I've fixed that. – Stéphane Laurent Jul 1 '17 at 15:24

Yeah, it's a multivariate t density.

Multiplying the three densities together and integrating (also using some properties of determinants) gives \begin{align*} &\iint \frac{|K|^{1/2}}{(2\pi)^{d/2}}e^{-0.5([X- \mu]^T K[X-\mu])} \frac{|\lambda_0K|^{1/2}}{(2\pi)^{d/2}}e^{-0.5([\mu - \mu_0]^T \lambda_0K[\mu-\mu_0])} \frac{|K| ^{(v-d-1)/2}}{2^{vd/2}|W|^{v/2}\Gamma_p|\frac{\nu}{2}|} e^{-tr(V^{-1}K)/2} d\mu dK \\ &= \iint\frac{|K|^{(v-d+1)/2}\lambda_0^{d/2} }{2^{d(1+\nu/2)} \pi^{d}\Gamma_p|\frac{\nu}{2}||W|^{v/2}} \exp\left[-\frac{Q}{2}\right] d\mu dK\\ \end{align*}

where \begin{align*} Q &= (X- \mu)^T K(X-\mu) + (\mu - \mu_0)^T \lambda_0K(\mu-\mu_0) + tr(V^{-1}K) \\ &= (\mu-X)^T K(\mu-X) + (\mu - \mu_0)^T \lambda_0K(\mu-\mu_0) + tr(V^{-1}K) \\ &= \mu^TK\mu + \mu^T\lambda_0K\mu- 2\mu^T(KX +\lambda_0K\mu_0)+ X^TKX + \mu_0^T\lambda_0K\mu_0 + tr(V^{-1}K) \\ &= (\mu - A)^TK(1 + \lambda_0)(\mu - A) - A^TK(1+\lambda_0)A+ X^TKX + \mu_0^T\lambda_0K\mu_0 + tr(V^{-1}K) \\ &= R - A^TK(1+\lambda_0)A+ X^TKX + \mu_0^T\lambda_0K\mu_0 + tr(V^{-1}K) \end{align*} with $A = \frac{K^{-1}}{(1 + \lambda_0)}(KX +\lambda_0K\mu_0) = (X + \lambda_0\mu_0)/(1+\lambda_0)$.

Completing the square like this helps us to recognize the normal density that integrates to $1$. After integrating out $\mu$, the last double integral simplifies to

$$\left(\frac{\lambda_0}{1+\lambda_0}\right)^{d/2} \int\frac{|K|^{(v-d)/2} }{2^{d(1+\nu)/2} \pi^{d/2}\Gamma_d|\frac{\nu}{2}||W|^{v/2} } \exp\left[-\frac{G}{2}\right] dK$$ where (using properties of the trace operator) \begin{align*} G &= - A^TK(1+\lambda_0)A+ X^TKX + \mu_0^T\lambda_0K\mu_0 + tr(V^{-1}K)\\ &= tr(-AA^TK(1+\lambda_0))+ tr(XX^TK) + tr(V^{-1}K) \\ &= -tr((X + \lambda_0\mu_0)(X + \lambda_0\mu_0)^TK)/(1+\lambda_0)^2+ tr(XX^TK) + tr(V^{-1}K) \\ &= tr\left( \left\{(X + \lambda_0\mu_0)(X + \lambda_0\mu_0)^T)/(1+\lambda_0)^2 + XX^T + V^{-1} \right\}K\right) \\ &= tr(C^{-1}K). \end{align*}

Now we can integrate out $K$ because we recognize the Wishart distribution \begin{align*} &\left(\frac{\lambda_0}{1+\lambda_0}\right)^{d/2}\pi^{-d/2}|W|^{-v/2} \frac{1}{\Gamma_d\left(\frac{\nu}{2}\right)} \int\frac{|K|^{([v+1]-d-1)/2} }{2^{d(\nu+1)/2} } \exp\left[-\frac{tr\left[C^{-1}K\right]}{2}\right] dK \\ &= \left(\frac{\lambda_0}{1+\lambda_0}\right)^{d/2}\pi^{-d/2}|W|^{-v/2} \frac{\Gamma_d\left(\frac{\nu+1}{2}\right)}{\Gamma_d\left(\frac{\nu}{2}\right)} |C|^{(\nu+1)/2}. \end{align*}

• Thank you very much Taylor, that's awesome. I really appreciate it. Did you do the derivation yourself or did you find it somewhere? – ejlouw Jul 1 '17 at 18:11
• @ejlouw I did this myself last night to refresh my memory on all the tricks, so I don't have any citations for you. I was looking mostly at Wiki pages. It's a pretty widely used result that's probably in a lot of Bayesian texts. See Stéphane Laurent's comment – Taylor Jul 1 '17 at 19:21

It is a well-known problem. A complete set of solutions can be found at https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf. If speed is of the essence then this is an acceptable solution, but there are concerns about using the Wishart prior. See, for example, https://dahtah.wordpress.com/2012/03/07/why-an-inverse-wishart-prior-may-not-be-such-a-good-idea/

edit note I updated the link to include a solution with a predictive distribution. Sorry, didn't notice.

• Just to elaborate on Dave's point (+1), the Gaussian-Wishart model is computationally convenient because you can exploit conjugacy in Gibbs sampling. But if you're not bottlenecked on computation, employing a prior over correlations and a prior over variance implies a model which allows each to vary independently. This strategy is feasible (and actually preferred for reasons of geometry) in HMC-NUTS. See the Stan documentation for more details – Sycorax Jun 30 '17 at 15:00
• Can we maybe see it typed up? You linked to two documents, and of them is nearly 30 pages. Also, just because a problem is well-known doesn't mean it's solveable. – Taylor Jun 30 '17 at 22:44
• @taylor the solution for the predictive density is in section 8 at equation 232. No, I won't type it out because much of the back work occurs early in the document and the amount of work isn't trivial to do. There is a reason this is a thirty-page document with appendices. – Dave Harris Jun 30 '17 at 23:11

The answer can be derived from the following result. If $\Sigma \sim {\cal IW}_\nu(V)$ (inverse-Wishart) and $(G \mid \Sigma) \sim {\cal N}(\theta, \lambda\Sigma)$, then $G \sim {\cal T}_{\nu-d+1}\left(\theta, \lambda\frac{V}{\nu-d+1}\right)$ (multivariate Student), where $d$ is the dimension.

In your problem: \begin{align} K^{-1} & \sim {\cal IW}(\nu, V), \\ (\mu \mid K) & \sim {\cal N}\bigl(\mu_0, {(\lambda_0 K)}^{-1} \bigr), \\ (X \mid \mu, K) & \sim {\cal N}\bigl(\mu, K^{-1} \bigr). \end{align}

Thus you have $$(X \mid K) \sim {\cal N}\left(\mu_0, \left(1+\frac{1}{\lambda_0}\right)K^{-1} \right)$$

Applying the result I mentioned, we get $$X \sim {\cal T}_{\nu-d+1}\left(\mu_0, \left(1+\frac{1}{\lambda_0}\right)\frac{V}{\nu-d+1} \right).$$

• Thanks Stéphane, where did you find this answer? Or did you derive it yourself? – ejlouw Jul 1 '17 at 18:17
• @ejlouw You can find this result in the pdf posted by Dave Harris, but it is a bit hidden by the context. By the way it should be given in the wiki page about the Normal-Wishart distribution, but currently it is not (the page only mentions that the marginal distribution is a multivariate Student, without giving the parameters). – Stéphane Laurent Jul 1 '17 at 18:23