# Marginal likelihood of bivariate Gaussian model

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.