## Programme to simulate from a mixture of two Gaussians
 
# Define the density for the two-component normal mixture
mix <- function(x, m1, m2, sd1, sd2, p){

  return(p*dnorm(x, mean=m1, sd=sd1) + (1-p)*dnorm(x, mean=m2, sd=sd2))
}

# set the parameters
p = 0.6
m1 = 0
m2 = 4
sd1=1
sd2=2

# plot the density
curve(mix(x, m1, m2, sd1, sd2, p), from=-5, 10)

n = 10000 # number of samples

u = runif(n)
mix.comp  = as.numeric(u < p)

sim = rep(NA, n)
for(i in 1:n){
  if(mix.comp[i]){
    sim[i] = rnorm(1, mean=m1, sd=sd1)
  } else {
    sim[i] = rnorm(1, mean=m2, sd=sd2)
  }
}

library(MASS)
truehist(sim)
curve(mix(x, m1, m2, sd1, sd2, p), from=-5, 10, add=T)

# alternative
numComp1 = sum(mix.comp)
numComp2 = n - numComp1

sim2 <- c(rnorm(numComp1, mean=m1, sd=sd1), rnorm(numComp2, mean=m2, sd=sd2))
truehist(sim2)
curve(mix(x, m1, m2, sd1, sd2, p), from=-5, 10, add=T)