question about joint and marginal distribution I have a question about how to compute the joint and marginal distribution.
Let $x$ and $y$ have the Gaussian densities.
$x \sim \mathcal{N}(m, P), \  y|x \sim \mathcal{N}(Hx, R)$
how to compute $p(x,y)$ and $p(y)$?
The answer is $y\sim \mathcal{N}(Hm, HPH^T + R)$, but why?
In fact, i see it from the Kalman filter.
Thanks very much!
 A: The answer to the first part of your question is
$$p(x,y)=p(y|x)p(x)=p(x|y)p(y)$$
And hence
$$p(y)=\int p(x,y)\mathrm dx=\int p(y|x)p(x)\mathrm dx$$
UPDATE.
Your assumptions differ a bit from the model that's usually assumed in Kalman filtering (so I guess the answer will not be in the form that you have mentioned).
The derivation of the marginal pdf is painful, so I will just outline the main idea and several steps.
So we have $x \sim \mathcal{N}(m, P), \  y|x \sim \mathcal{N}(Hx, R)$ and $p(x,y)=p(y|x)p(x)$
Then  
$$
\begin{equation}
p(x,y)\propto \exp\left(-{(x-m)^TP^{-1}(x-m)}\right)\exp\left({-(y-Hx)^TR^{-1}(y-Hx)}\right) \qquad \text{(1)}
\end{equation}
$$
hereafter the incidental terms (that do not depend on $x$ or $y$) will be omitted.
 To get the marginal pdf $p(y)$ one needs to integrate $(1)$ with the respet to $x$, so one needs to complete the full squares in the exponents.
 Combining the exponent terms:
$$
p(x,y)\propto \exp{\left[-\left((x-m)^TP^{-1}(x-m)+(y-Hx)^TR^{-1}(y-Hx)\right)\right]} 
$$
Now I will work only with the power (without the minus sign):
$$x^TP^{-1}x-x^TP^{-1}m-m^TP^{-1}x+m^TP^{-1}m+y^TR^{-1}y-(Hx)^TR{-1}y-y^TR^{-1}(Hx)+(Hx)^TR^{-1}(Hx)$$
Which is equal to 
$$
x^T\left(P^{-1}+H^TR^{-1}H\right)x-2x^T\left(P^{-1}m+H^TR^{-1}y\right)+y^TR^{-1}y
+m^TP^{-1}m \qquad \text{(2)}$$
Here I used the fact that  $R$ are symmetric and so $(R^{-1})^T=R^{-1}$, the same is with $P$. And that $$(Hx)^T=x^TH^T$$ and that $$x^TP^{-1}x+x^TH^TR^{-1}Hx=x^T\left(P^{-1}+H^TR^{-1}H\right)x$$
Again I will omit the last term $m^TP^{-1}m$.
So we would like to complete the full square like:
$$(x-\mu_x)^T\Sigma_x(x-\mu_x)=x^T\Sigma_xx-2x^T\Sigma_x\mu_x+\mu_x^T\Sigma_x\mu_x$$
Comparing the last equation with $(2)$ one can find that:
$$\Sigma_x=\left(P^{-1}+H^TR^{-1}H\right)$$
$$\mu_x=\Sigma_x^{-1}\left(P^{-1}m+H^TR^{-1}y\right)$$
Adding and substracting $\mu_x^T\Sigma_x\mu_x$ (up to the omitted term $m^TP^{-1}m$) $(2)$ will look like:
$$(x-\mu_x)^T\Sigma_x(x-\mu_x)-\mu_x^T\Sigma_x\mu_x+y^TR^{-1}y$$
Now we have separated the part containing $x$ and hence we can integrate with the respect to $x$. Using the fact that $\int \exp{(-(x-\mu_x)^T\Sigma_x(x-\mu_x))} \mathrm dx$ equals to normalizing constant one can obtain the marginal pdf of $y$:
$$p(y)\propto \exp{(-y^TR^{-1}y+\mu_x^T\Sigma_x\mu_x)}$$
So it is gaussian, but to find its mean and variance one needs to perform some even more painful derivations.
A: First, you show that $(X, Y)$ is jointly Gaussian. Suppose $X$ is $k$-dimensional and $Y$ is $\ell$-dimensional.
\begin{eqnarray*}
p_{X,Y}(x,y) & = & p_X(x)p_{Y|X}(y,x) \\
& = & (2\pi)^{-k/2}|P|^{-1/2}\exp\left(-\frac 12(x - m)^T P^{-1} (x - m)\right) \cdot \\
& & (2\pi)^{-\ell/2}|R|^{-1/2}\exp\left(-\frac 12(y - Hx)^T R^{-1} (y - Hx)\right) \\
& = & (2\pi)^{-(k+\ell)/2}|PR|^{-1/2}\exp\left(-\frac 12(x - m)^T P^{-1} (x - m)\right. \\
& & \left.-\frac 12(y - Hx)^T R^{-1} (y - Hx)\right).
\end{eqnarray*}
We will be done if we can write the expression inside $\exp$ in the form
$$
-\frac 12\left(\begin{bmatrix}x^T & y^T\end{bmatrix} - \mu^T\right)\Sigma^{-1}
\left(\begin{bmatrix}x \\ y\end{bmatrix} - \mu\right)
$$
for some positive definite matrix $\Sigma^{-1}$ and vector $\mu$. This can be done easily, but the process is rather tedious. If I may suggest, start with writing $y - Hx$ as $(y - Hm) - H(x - m)$. You will find that
$$
\mu = \begin{bmatrix}
m \\ Hm
\end{bmatrix}
$$
and
$$
\Sigma^{-1} =
\begin{bmatrix}
P^{-1} + H^TR^{-1}H & H^TR^{-1} \\
R^{-1}H & R^{-1}
\end{bmatrix}
$$
work.
The non-trivial part is inverting $\Sigma$. You can do this directly by using Schur complement, or you can compute $Cov(X, Y)$ and $Var(Y, Y)$ from their definitions and use the fact that $\Sigma$ is made from these covariance matrices. I will demonstrate the second method here:
$$
E[Y] = E[E[Y\mid X]] = E[HX] = HE[X] = Hm.
$$
\begin{align*}
Cov(X, Y) & = E\left[(X - E[X])(Y - E[Y])^T\right] \\
& = E\left[E\left[(X - E[X])(Y - E[Y])^T\mid X\right]\right] \\
& = E\left[E\left[(X - m)(Y - Hm)^T\mid X\right]\right] \\
& = E\left[(X - m)(HX - Hm)^T\right] \\
& = E\left[(X - m)(X - m)^T\right]H^T \\
& = PH^T \\
\end{align*}
\begin{align*}
Var(Y) & = E\left[(Y - E[Y])(Y - E[Y])^T\right] \\
& = E\left[E\left[(Y - E[Y])(Y - E[Y])^T\mid X\right]\right] \\
& = E\left[E\left[((Y - HX) - H(m - X))((Y - HX) - H(m - X))^T\mid X\right]\right] \\
& = E\left[R + H(m - X)(m - X)^TH^T\right] \\
& = R + HPH^T.
\end{align*}
Therefore, $(X, Y)$ is jointly Gaussian with mean $\mu = \begin{bmatrix}m\\Hm\end{bmatrix}$ and variance $\Sigma = \begin{bmatrix}P & PH^T \\ HP & R + HPH^T\end{bmatrix}$.
To find the marginal distribution of $Y$, just take the submatrix of $\Sigma$ that corresponds to $Y$.
