Suppose $$x=Pt+e$$ $$y=Ct+f$$ where $t\sim N(0,I)$, $e\sim N(0,\Sigma_x)$, and $f\sim N(0,\Sigma_y)$. Also, $\Sigma_x\in R^{d\times d}$, $\Sigma_y\in R^{k\times k}$, and $t\in R^r$ One can compute the marginal probability density function of $(x,y)$ as following $$p(x,y)=\int p(x,y,t)dt=\int N(Pt,\Sigma_x)N(Ct,\Sigma_y)N(0,I)dt$$ Expanding out and rearranging the terms in the exponents, the above equation becomes the following $$\dfrac{1}{2\pi^{(d+k+r)/2}|\Sigma_x|^{1/2}|\Sigma_y|^{1/2}}\int \exp(x^T\Sigma_x^{-1}x+y^T\Sigma_y^{-1}y-2(x^T\Sigma_x^{-1}P^T+y^T\Sigma_y^{-1}C^T)t+t^T(P^T\Sigma_x^{-1}P+C^T\Sigma_y^{-1}C+I)t)dt$$ I realize that the complex term in the exponent of the above equation is a quadratic form of $t$ which might be rearranged into the form $(t-\mu_t)^T\Sigma_t^{-1}(t-\mu_t)$. I managed to obtain $\Sigma_t=(P^T\Sigma_x^{-1}P+C^T\Sigma_y^{-1}C+I)^{-1}$. Is it therefore correct to simplify the integral as the following by using the property of the multivariate Gaussian distribution? $$\int \exp(x^T\Sigma_x^{-1}x+y^T\Sigma_y^{-1}y-2(x^T\Sigma_x^{-1}P^T+y^T\Sigma_y^{-1}C^T)t+t^T(P^T\Sigma_x^{-1}P+C^T\Sigma_y^{-1}C+I)t)dt=2\pi^{r/2}|P^T\Sigma_x^{-1}P+C^T\Sigma_y^{-1}C+I|^{1/2}$$
If it is correct, the marginal $p(x,y)$ becomes irrelevant of $(x,y)$ which is weird.Can anyone point out where I might got it wrong?
Update and Correction $$$$ I have now identified the problem with the above derivation, which was all due to a fundamental mistake as explained below. I have assumed, due to $p(x,y,t)\propto p(t|x,y)$ $$p(x,y)=\int p(x,y,t)dt=\int p(t|x,y)dt$$ which is absurdly wrong. The correct way of doing it is through partitioning. $$P(x,y,t)=\begin{bmatrix}\begin{bmatrix}x^T\\y^T\end{bmatrix}\\t^T\end{bmatrix}\begin{bmatrix}\begin{bmatrix}\Sigma_x^{-1}&0\\0&\Sigma_y^{-1}\end{bmatrix}&\begin{bmatrix}\Sigma_x^{-1}P^T\\\Sigma_y^{-1}C^T\end{bmatrix}\\\begin{bmatrix}P\Sigma_x^{-1}&C\Sigma_y^{-1}\end{bmatrix}&P^T\Sigma_x^{-1}P+C^T\Sigma_y^{-1}C\end{bmatrix}\begin{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}\\t\end{bmatrix}$$
One can easily compute $Var[x,y]$ by applying the matrix inversion lemma on the square block matrix in the middle of the above expression. It is also easy to see that $E(x,y)=0$.