Sampling from Gaussian mixture models, when are the sampled data independent? Suppose I generate a Gaussian mixture model with $N$ Gaussian distributions
$p(x) = \sum\limits_{i = 1}^N w_i \mathcal{N}(x;\mu_i, \Sigma_i)$
where $w_i$ are the weights.
Now I sample some points $\{x_n\}$ from $p(x)$
What condition is needed to ensure that the points are independent from each other?
Is it sufficient to assume that each component $\mathcal{N}(x;\mu_i, \Sigma_i)$ have diagonal covariance matrix $\Sigma_i$?
 A: 
Main point: You are confusing independence between the simulated
  vectors and independence between the components of a single vector.

Unless otherwise stated, two simulations $X_1,X_2$ from$$\sum\limits_{i = 1}^N w_i \mathcal{N}(x;\mu_i, \Sigma_i)\tag{1}$$will be independent, while the components of the vector $X_1$ will most likely be dependent, even when the matrices $\Sigma_i$ are diagonal. (The intuition about that last point is that if $X_{11}$ is closer to $\mu_{1k}$ than to the other Normal means, then the other components of $X_1$ are also more likely to originate from this same $k$-th element of the mixture.) For instance a most standard way to simulate from (1) is to simulate
$$Z_1\sim\mathcal M(1;w_1,\ldots,w_N)\qquad X_1|Z_1=k \sim \mathcal{N}(x;\mu_k, \Sigma_k)\tag{2}$$and$$Z_2\sim\mathcal M(1;w_1,\ldots,w_N)\qquad X_2|Z_2=k \sim \mathcal{N}(x;\mu_k, \Sigma_k)$$If the four variables are mutually independent then $X_1$ and $X_2$ are independent.
As detailed in this other answer, the covariance matrix of the mixture is given by
$$\sum_{i=1}^N w_i \Sigma_i + \sum_{i=1}^N w_i(\mu_i-\bar\mu)(\mu_i-\bar\mu)^T$$when $\bar\mu=w_1\mu_1+\cdots w_N\mu_N$. This explains why the means play a role in the correlation. Another explanation is the above latent variable decomposition (2): conditionally on $Z_1=k$, the components of $X_1$ are independent if $\Sigma_k$ is diagonal, but the independence does not stay unconditionally, because all components of $X_1$ share the same $Z_1$.
Here is a short experiment to demonstrate the dependence: I simulated 10³ realisations from
$$\frac{1}{3}\mathcal N_2\left(\left[\matrix{1\\ 2\\} \right], \left[\matrix{1 &0\\0 &2\\}\right]\right)+\frac{1}{3}\mathcal N_2\left(\left[\matrix{1\\ 0\\} \right], \left[\matrix{2 &0\\0 &1\\}\right]\right)+\frac{1}{3}\mathcal N_2\left(\left[\matrix{-2\\ -2\\} \right], \left[\matrix{1 &0\\0 &1\\}\right]\right)$$
and obtained the following sample
$\qquad\qquad\qquad$
with a regression line over $x_2$ against $x_2$ corresponding to a correlation of $0.46$, despite the three matrices $\Sigma_i$ being diagonal.
