I assume the following model for a sample $y_1 \in \mathbb{R}^2$ of size $1$ with bivariate Gaussian likelihood and independent bivariate Gaussian and inverse-Wishart prior for the mean and variance parameters, respectively. That is
$$ y_1 \mid \mu, \Sigma \sim \mbox{Normal}(y_1 \mid\texttt{mean}=\mu, \texttt{Var}=\Sigma) \\ \mu \sim \mbox{Normal}(\mu \mid\texttt{mean}=\mu_\mu, \texttt{Var}=\Sigma_\mu)\\ \Sigma \sim \mbox{Inverse-Wishart}(\Sigma\mid\texttt{df}=\nu_\Sigma,\texttt{scale}=S_\Sigma) $$
I want to derive the analytical expression of the marginal likelihood of $y_1$, that is $$ p(y_1)=\iint p(y_1 \mid \mu, \Sigma) \, p(\mu) \, p(\Sigma) \mathrm{d}\mu \mathrm{d}\Sigma, $$ where, with an abuse of notation, $p(x)$ denotes the density of $x$ evaluated in the point $x$.
Remark: I know that if I assume a conjugate Normal-Inverse-Wishart prior such result is well-known and coincide with a Generalized-T distribution evaluated in $y_1$. See for instance Wikipedia Conjugate Prior. However, I do not want to change my prior.
My attempt: Relying on well-known conjugacy results (see Wikipedia Conjugate Prior) I know how to solve one of the two integral fixing either $\mu$ or $\Sigma$, but I cannot solve (that is find a tractable analytical expression in terms of well-known functions) the second integral.
