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
 A: 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.
A: 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*}
A: 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).
$$
