I would like to understand how to integrate out one of the dimensions of a multi-dimensional Gaussian. Here is a simple setting:
x ~ N($\mu_x$,$\sigma_x^2$) y ~ N($\mu_y$,$\sigma_y^2$)
We form a new random variable through linear addition: z = ax +by
Then z ~ N($a\mu_x + b\mu_y$,$a^2\sigma_x^2 + b^2\sigma_y^2$)
Now suppose I want to get back the marginal pdf of x. Then: $f_x(x) = \int_{-\infty}^{\infty} f(x,y)dy$.
More explicitly, $f_x(x) = \frac{1}{\sqrt{2\pi(a^2\sigma_x^2 + b^2\sigma_y^2)}}\int_{-\infty}^{\infty} e^{\frac{(ax+by-a\mu_x + b\mu_y)^2}{2(a^2\sigma_x^2 + b^2\sigma_y^2)}}dy$.
One option is to expand the exponent. Some terms come out of the integral, but the remaining terms contain a number of cross terms remaining. Specifically: 2abxy -2abx$\mu_x$ - 2aby$\mu_y$ +2ab$\mu_x\mu_y$.
How is this normally handled?
If I understand your question, you have a multidimensional gaussian rv $X=(x_1,\dots, x_n)$, and you would like to find the distribution of $(x_2,\dots, x_n)$. In general, if $X$ is multidimensional gaussian with mean $\mu$ and covariance $C$, and $A: \mathbb{R}^n\to\mathbb{R}^m$ is any linear map, then $AX$ is also multidimensional gaussian with mean $A\mu$ and covariance $A CA^T$ (this can be seen most easily using the formula for the characteristic function of a gaussian). In your case, the matrix $C$ is diagonal with entries $\sigma_x^2,\sigma_y^2$ and $A$ is the matrix of the projection onto the first coordinate $A=(1,0)$, thus recovering the distribution of X