Distribution of the sum of the two dependent bivariate gaussian distributions and related questions This is something I was thinking about and I decided to modify a question from a mid-term to ask this. 
Suppose $X_{1}$ and $X_{2}$ are two bivariate gaussian variables, decribed as 
$$
X_{i}=\begin{bmatrix}
        \mu_{1}+y\\
        \mu_{2}+w\\
        \end{bmatrix}
        $$
where $y$ and $w$ are iid random variables each with mean $0$ and variance $1$. $\mu_{1}$ and $\mu_{2}$ are unknown parameters.
The pdf of $X_{i}$ is $N(\mu,\Sigma)$, where $$\mu=\begin{bmatrix}
        \mu_{1}\\
        \mu_{2}\\
        \end{bmatrix}
        $$
and $$\Sigma=\begin{bmatrix}
        1&0\\
        0&1\\
        \end{bmatrix}
        $$
So, finding the joint pdf of $X_{1}$ and $X_{2}$, if they are independent, is easy enough, and may be generalised to n iid variables.
Firstly, what would the pdf be $X_{1}$ if y and w weren't iid? I should simply replace the off-diagonal entries of  $\Sigma$ with the co-variance, correct?
Now, suppose $X_{1}$ and $X_{2}$ aren't iid. What would the pdf of $Z=X_{1}+X_{2}$ be (assuming that this sum is also gaussian)? I have read around on stackexchange and there is an answer to this question if $X_{1}$ and $X_{2}$ were univariate, but I don't necessarily see how the maths would work out for my case.
Now, I also want to ask how I can put this in the context of a machine learning problem. Does $X_{i}$ being a vector mean that each member of the vector is a feature of the sample $X_{i}$? When are features of a sample dependent?
What would dependent $X_{i}$ physically mean, then? Can I model a system which collects samples as this: "If the first sample $X_{1}$, modelled as a Gaussian random variable, is less than 5000, then $X_{2}$ is obtained iid. Otherwise, the mean of $X_{2}$ shifts by an amount $a$."  
I hope that this question isn't too vague or worded confusingly.  
 A: If you ask "What if they're not independent?" that doesn't say what they are instead.  There are numerous possibilities. You ask whether you should "simply replace the off-diagonal entries of $\Sigma$ with the co-variance". The answer is yes if the pair $(x,y)$ is jointly normally distributed.  What that means is that it is so distributed that every linear combination $ax+by$ is normally distributed.  But there are joint distributions in which $(x,y)$ is not jointly normal and yet $x$ and $y$ are separately normal.  In some such cases $x$ and $y$ are correlated, and there is a "covariance matrix" in which instead of $0$ one has the covariance between the two, but that doesn't specify the distribution --- more information is needed after that.  And in the case where they are jointly normal, the mean vector and covariance matrix don't specify the distribution without the additional information that they are jointly normal.
So: If you have the information saying they are jointly normal rather than only separately normal, then specifying the mean vector and covariance matrix is enough.  To say that more precisely, suppose $(X,Y)$ is a jointly normally distributed pair of random variables and $(U,V)$ is also a jointly normally distributed pair of random variables, and both pairs have the same mean vector and the same covariance matrix.  Then the distribution of $(X,Y)$ is the same as the distribution of $(U,V)$.
You then ask what if $X_1, X_2 \in\mathbb R^2$ are not independent.  Again we have the question of what they are instead.  It might be that $(X_1,X_2)\in\mathbb R^4$ is a jointly normally distributed quadruple of real random variables.  Alternatively, it might not, and then a lot more information is needed to specify how it is distributed.  But let us suppose
$$
\begin{bmatrix} X_1 \\ X_2 \end{bmatrix} \sim N(\mu,\Sigma) \text{ and } \Sigma = \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{12}^\top & \Sigma_{22} \end{bmatrix}
$$
where each $\Sigma_{ij}$ is a $2\times2$ real matrix.  Then we have $\operatorname{E}(X_1+X_2) = \operatorname{E}(X_1)+\operatorname{E}(X_2)$. Below I will use $\text{“}{\operatorname{var}}\text{''}$ to mean the "covariance matrix" (following William Feller's famous book on probability).
$$
\operatorname{var}(X_1+X_2) = \Sigma_{11} + \Sigma_{22} + \Sigma_{12} + \Sigma_{12}^\top,
$$
and the sum $X_1+X_2\in\mathbb R^2$ is jointly normally distributed (in the sense that every linear combination of its two components has a univariate normal distribution), and there is only one bivariate normal distribution that has that particular mean vector and that particular covariance matrix.
The proofs of the assertions above take more work that what is now in this posted answer.
