Find the variance-covariance matrix for a linear combination of multiple bivariate normal distribution? I have an arbitrary number of independnet bivariate normal distributions with $\mu_i = [x_i,z_i]$ & $ \Sigma_i=
\left(\begin{array}{cc} 
\sigma^2_{x_i} & \sigma^2_{x_i,z_i}\\\ 
\sigma^2_{x_i, z_i} & \sigma^2_{z_i}
\end{array}\right)
$
Where i is arbitrarily large
I want to take a linear combination of these bivariate normal distributions with weights $c = [c_1,...,c_i]$ where $\sum c_i = 1$ & $c_i >0$
Obviously, the linear combination of $\mu_{mixture} = [\sum c_ix_i,\sum c_iz_i]$
However, I am not sure about the linear combination of the variance-covariance matrix.
Does anyone know how I can calculate this pooled & weighted variance-covariance? Looking for the variance-covariance matrix for the mixture distribution.
Thanks so much!

Let $X_i\stackrel{\text{indep}}\sim\mathcal N(\mu_i,\Omega_i)$ and let $S = \sum_{i=1}^n c_iX_i$. A linear combination of independent Gaussians is Gaussian so we just need the mean and covariance. By linearity we have
$$
\text E[S] = \sum_i c_i\mu_i
$$
and by independence we have
$$
\text{Var}[S] = \sum_i \text{Var}[c_iX_i] = \sum_i c_i^2 \Omega_i
$$
so
$$
S\sim\mathcal N\left(\sum_i c_i\mu_i, \sum_i c_i^2\Omega_i\right).
$$
This applies no matter what the $c_i$ are and for any dimension of $X_i$.

The above part assumed $n < \infty$. If we have a countably infinite number of $X_i$ then whether or not the series $\sum_{i=1}^\infty c_i X_i$ converges depends on how the $c_i$, $\mu_i$, and $\Omega_i$ evolve and we can use Kolmogrov's three series theorem to understand when this happens.

I interpreted this to mean you wanted the distribution of a linear combination of Gaussians. If you meant a finite mixture of Gaussians then we can work it out in the following way. Let $f_i$ be the density of $X_i$ and let $S \sim \sum_{i=1}^n c_i f_i$ be the mixture distribution. You didn't state that $c_i \geq 0$ but I'll assume that so that this is a valid density. Then we have
$$
\text E[S] = \int s \sum_i c_i f_i(s)\,\text ds = \sum_i c_i \text E[X_i] = \sum_i c_i \mu_i
$$
as before, except now this represents a convex combination of the $\mu_i$ where that was not guaranteed before. I'll use $\mu_\text{mix} = \sum_i c_i\mu_i$ as the mixture mean.
For the variances we need
$$
\text E[SS^T] = \int ss^T \sum_i c_i f_i(s)\,\text ds = \sum_i c_i \text E[X_iX_i^T]
$$
so all together
$$
\text{Var}[S] = \text E[SS^T] - (\text E S)(\text ES)^T \\ =\sum_i c_i \text E[X_iX_i^T] - \mu_\text{mix}\mu_\text{mix}^T
$$
which is more complicated than $\sum_i c_i^2\Omega_i$
