Simulate from a truncated mixture normal distribution I want to simulate a sample from a mixture normal distribution such that
$$p\times\mathcal{N}(\mu_1,\sigma_1^2) + (1-p)\times\mathcal{N}(\mu_2,\sigma_2^2) $$
is restricted to the interval $[0,1]$ instead of $\mathbb{R}$. This means I want to simulate a truncated mixture of normal distributions. 
I know that there are some algorithm to simulate a truncated normal (i.e from this question) and corresponding package in R to do this. But how can I simulate a truncated mixture normal? Is it the same if I simulate two truncated normal of $\mathcal{N}(\mu_1,\sigma_1^2)$ and $\mathcal{N}(\mu_2,\sigma_2^2$) to make a truncated mixture normal? 
 A: Simulation from a truncated normal is easily done if you have access to a proper normal quantile function. For instance, in R, simulating $$
\mathcal{N}_a^b(\mu,\sigma^2)$$where $a$ and $b$ denote the lower and upper bounds can be done by inverting the cdf $$\dfrac{\Phi(\sigma^{-1}\{x-\mu\})-\Phi(\sigma^{-1}\{a-\mu\})}{\Phi(\sigma^{-1}\{b-\mu\})-\Phi(\sigma^{-1}\{a-\mu\})}
$$
e.g., in R
x = mu + sigma * qnorm( pnorm(a,mu,sigma) + 
     runif(1)*(pnorm(b,mu,sigma) - pnorm(a,mu,sigma)) )

Otherwise, I developed a truncated normal accept-reject algorithm twenty years ago.
If we consider the truncated mixture problem, with density
$$
f(x;\theta) \propto \left\{p\varphi(x;\mu_1,\sigma_1)+(1-p)\varphi(x;\mu_2,\sigma_2)\right\}\mathbb{I}_{[a,b]}(x)
$$
it is a mixture of truncated normal distributions but with different weights:
$$
f(x;\theta) \propto p\left\{\Phi(\sigma_1^{-1}\{b-\mu_1\})-\Phi(\sigma_1^{-1}\{a-\mu_1\}) \right\}\dfrac{\sigma_1^{-1}\phi(\sigma_1^{-1}\{x-\mu_1\})}{\Phi(\sigma_1^{-1}\{b-\mu_1\})-\Phi(\sigma_1^{-1}\{a-\mu_1\})}  \\[15pt]
+(1-p)\left\{\Phi(\sigma_2^{-1}\{b-\mu_2\})-\Phi(\sigma_2^{-1}\{a-\mu_2\}) \right\}\dfrac{\sigma_2^{-1}\phi(\sigma_2^{-1}\{x-\mu_2\})}{\Phi(\sigma_2^{-1}\{b-\mu_2\})-\Phi(\sigma_1^{-1}\{a-\mu_2\})}
$$
Therefore, to simulate from a truncated normal mixture, it is sufficient to take
$$x=\begin{cases}
x_1\sim\mathcal{N}_a^b(\mu_1,\sigma_1^2) &\text{with probability }\\
&\qquad p\left\{\Phi(\sigma_1^{-1}\{b-\mu_1\})-\Phi(\sigma_1^{-1}\{a-\mu_1\}) \right\}\big/\mathfrak{s}\\
x_2\sim\mathcal{N}_a^b(\mu_2,\sigma_2^2) &\text{with probability }\\
&\qquad(1-p)\left\{\Phi(\sigma_2^{-1}\{b-\mu_2\})-\Phi(\sigma_2^{-1}\{a-\mu_2\}) \right\}\big/\mathfrak{s}
\end{cases}
$$
where
\begin{align}
\mathfrak{s}=&p\left\{\Phi(\sigma_1^{-1}\{b-\mu_1\})-\Phi(\sigma_1^{-1}\{a-\mu_1\}) \right\}+  \\
&(1-p)\left\{\Phi(\sigma_2^{-1}\{b-\mu_2\})-\Phi(\sigma_2^{-1}\{a-\mu_2\}) \right\}
\end{align}
