Q:

Two Correct Ways to Sample from a Two Component Mixture but Different Variance Results...What's happening?

Let $X_1 \sim \operatorname{Poisson}(7)$ and $X_2 \sim \operatorname{Poisson}(3)$ be two independent random variables. 
Let $\{0.7,0.3\}$ be our collection of mixing weights.
The following R code is to sample from the mixture distribution of X1 and X2
    mixture_poisson <- function(n,lambdas,weights) 
    {
     I = sample(1:length(weights),size=n,replace=T,prob=weights)
     mixture = rpois(n,lambda[I])

     mean(mixture)
     var(mixture)
    }

The problem that I am having is that the variance from the experiment is much higher than the theoretical variance, which should be 
$$7(0.7)^2 + 3(0.3)^2 = 3.7$$
Every time I run the program, I get variance outputs like $9.42$ or $9.41$. 
However, if I change the code to something like
    mixture_poisson <- function(n,lambdas,weights) 
    {
     mixture = weights[1]*rpois(n,lambdas[1]) + weights[2]*rpois(n,lambdas[2])

     mean(mixture)
     var(mixture)
    }

Then I have great approximations for the variance, with values like $3.85$ or $3.59$.
So my question is, why does the first block of code return horrendous estimates of the variance? What is the first block of code doing that's different from the second block of code, because to my mind they are doing the same thing.
Any insight would be great! Thank you!!

A:

as it happens, both proposals in the question are incorrect:

in the first proposal the simulation is correctly implemented by sample but [there is a bug in the R code and] the computation of the variance of a mixture $$f(x)=p_1 f_1(x)+ p_2 f_2(x)$$is incorrect: it should be $$p_1 \text{var}_{f_1}(X)+p_2\text{var}_{f_2}(X)+p_1\mathbb{E}_{f_1}(X)^2+p_2\mathbb{E}_{f_2}(X)^2-\mathbb{E}_{f}(X)^2$$which leads to
$$p_1\lambda_1+p_2\lambda_2+p_1\lambda_1^2+p_2\lambda_2^2-(p_1\lambda_1+p_2\lambda_2)^2$$which gives a numerical value of $9.16$ for the reference parameters
in the second proposal the simulation is invalid as it returns a weighted average of two Poisson samples, which is not distributed from a Poisson mixture but from a Poisson convolution, in which case the variance is $$p_1^2\lambda_1+p_2^2\lambda_2$$

