How to merge two Gaussians I have two multivariate Gaussians each defined by mean vectors and Covariance matrices (diagonal matrices). I want to merge them to have a single Gaussian i.e. I assume there is only one Gaussian but I separated observations randomly into two groups to get two different Gaussians which are not too different than each other.
Since I know the number of observations in each of two Gaussians, combined mean estimation is straight forward : $\frac{n_1\mu_1 + n_2\mu_2}{n_1+n_2}$
But, what about the Covariance matrix?
Thanks   
EDIT:
The question was confusing in the original post, especially the "merging Gaussians" part. Maybe the following paragraph would be a better choice.
I have two sets of observations drawn from two multivariate Gaussians each defined by mean vectors and Covariance matrices (diagonal matrices). I want to merge the observations to have a single sample, and I assume to have another Gaussian (i.e. I assume initially there was only a single Gaussian, and observations were separated into two groups to get two different Gaussians).
 A: Given two multivariate normal distributions $N(\mu_0, \Sigma_0)$ and $N(\mu_1, \Sigma_1)$, the distribution you want is $N(\mu_f,\Sigma_f)$,
where
$\Sigma_f = \Sigma_1 (\Sigma_1 + \Sigma_0)^{-1}\Sigma_0$
$\mu_f = \Sigma_1 (\Sigma_1 + \Sigma_0)^{-1} \mu_0 + \Sigma_0 (\Sigma_1 + \Sigma_0)^{-1} \mu_1$
This works for any covariance matrices (no need to be diagonal).
A: Ok I solved it :)
Since covariance matrix is diagonal we can assume having multiple univariates. And then variance combination is as
$$\hat{\mu} = \frac{n_1\mu_1 + n_2\mu_2}{n_1+n_2}$$
$$\hat{\sigma}^2 = \frac{(\sigma_1^2 + \mu_1^2)n_1 + (\sigma_2^2 + \mu_2^2)n_2}{ (n_1+n_2)} - \hat{\mu}^2$$
Here, I used $\sigma^2 = E[x^2] - E[x]^2$
thanks again
A: As the answer of JulienD is good. But I would like to rewrite the final formulas in a better way and extend them to many random variables, and want to clarify that the approved answer is likely not correct.
First, I rewrite the answer of JulienD with few additional steps to give it a better form as follows:

*

*Single variate:
\begin{align}
\mu &= E[x] = \frac{n_1\mu_1 + n_2\mu_2}{n_1 + n_2}
\end{align}
\begin{align}
\sigma^2 & = E[x^2] - E[x]^2 \\
& = \frac{\sigma_1^2 n_1^2 + \sigma_2^2 n_2^2 + (n_1\mu_1 + n_2\mu_2)^2}{(n_1+n_2)^2} - \mu^2 \\
& = \frac{\sigma_1^2 n_1^2 + \sigma_2^2 n_2^2}{(n_1+n_2)^2} + \frac{(n_1\mu_1 + n_2\mu_2)^2}{(n_1+n_2)^2} - \mu^2 \\ 
& = \frac{\sigma_1^2 n_1^2 + \sigma_2^2 n_2^2}{(n_1+n_2)^2} + \mu^2 - \mu^2 \\
& = \frac{n_1^2}{(n_1+n_2)^2}\sigma_1^2  + \frac{n_2^2}{(n_1+n_2)^2}\sigma_2^2 
\end{align}

*

*Multivariate:
\begin{align}
\Sigma &= E[XX^T] - E[X]E[X]^T \\
& = \frac{n_1^2 \Sigma_1 + n_2^2 \Sigma_2 + (n_1\mu_1 + n_2\mu_2)(n_1\mu_1 + n_2\mu_2)^T}{(n_1 + n_2)^2} - \mu\mu^T \\
&= \frac{n_1^2 \Sigma_1 + n_2^2 \Sigma_2}{(n_1 + n_2)^2} + \frac{(n_1\mu_1 + n_2\mu_2)(n_1\mu_1 + n_2\mu_2)^T}{(n_1 + n_2)^2} - \mu\mu^T \\
& = \frac{n_1^2 \Sigma_1 + n_2^2 \Sigma_2}{(n_1 + n_2)^2} + \mu\mu^T - \mu\mu^T \\
&= \frac{n_1^2}{(n_1 + n_2)^2} \Sigma_1 + \frac{n_2^2}{(n_1 + n_2)^2} \Sigma_2
\end{align}
Another approach to get a formula of variance or covariance matrix (for multivariate) is that we use the theorem of variance (which can find in many textbooks of statistics, e.g., Larry Wasserman, All of Statistics, theorem 3.20, page 52, or prove via definition as a similar way as JulienD did):
$$
\mathbb{V}(X+Y) = \mathbb{V}(X) + \mathbb{V}(Y) + 2 \textbf{Cov}(X,Y)
$$
More generally, for random variables $X_1, X_2,..., X_N$:
$$
\mathbb{V}(\sum_{i=1}^{N}a_i X_i) = \sum_{i=1}^{N}a_i^2\mathbb{V}(X_i) + 2 \sum_{i=1}^{N-1}\sum_{j=i+1}^{N} a_i a_j \textbf{Cov}(X_i,X_j)
$$
Consider $X_i \sim \mathcal{N}(\mu_i, \Sigma_i), i=1 \cdots N,$; and they are independent of each other. Hence, the covariances are zero, i.e., $\textbf{Cov}(X_i,X_j)=0.$
A newly aggregated these Gaussian distributions is defined as a weighted sum:
$$X=\sum_{i=1}^{N} a_i X_i=\sum_{i=1}^{N} \frac{n_i}{\sum_{l=1}^{N}n_l} X_i,$$
where $a_i=\frac{n_i}{\sum_{l=1}^{N}n_l}$ and $\sum_{i=1}^N a_i = 1$.
So, we get the mean and covariance of the aggregated distribution as follows:
$$
\mu = \sum_{i=1}^{N} a_i \mu_i
$$
\begin{align}
\Sigma &=\mathbb{V}(X)= \mathbb{V}(\sum_{i=1}^{N} a_i X_i)\\
&=\sum_{i=1}^{N}a_i^2\mathbb{V}(X_i) + 2 \sum_{i=1}^{N-1}\sum_{j=i+1}^{N} a_i a_j \textbf{Cov}(X_i,X_j) \\
& = \sum_{i=1}^{N}a_i^2 \Sigma_i + 0 \\
& = \sum_{i=1}^{N}a_i^2 \Sigma_i
\end{align}
Now, look at the answer of ahmethungari, it seems the answer is not correct.
