pdf1 = function(n)
{
# Generate n i.i.d. random variates from
#
# f(x) = 0.2 N1(x) + 0.3 N2(x) + 0.5 N3(x)
#
# N1: normal(0   , 0.50)
# N2: normal(6.5 , 0.25)
# N3: normal(14.5, 0.75)
#
# version 1: not the most memory-efficient!

  # 1. generate the indexes (1,2,3) of the 
  #    mixture components, with probabilities
  #    0.2, 0.3 and 0.5, respectively:
  idx = rep(2, n)
  u   = runif(n) 
  idx[u < 0.2] = 1
  idx[u >= 0.5] = 3

  # 2. generate samples from the 3 normal populations
  mxt = rbind(
      rnorm(n, mean=0, sd=0.5),
      rnorm(n, mean=6.5, sd=1.25),
      rnorm(n, mean=14.5, sd=0.75)
      )
    
  # 3. select the component of the mixture according
  #    to priors in idx
  i = idx + (seq(1,1000)-1)*3
  x = mxt[i]

  return(x)
}

# x = pdf1(1000)
# plot(density(x), ylim=c(0,0.15))
# xx = seq(-6,22,.1)
# y = .2*dnorm(xx,0,0.5) + .3*dnorm(xx,6.5,1.25) + .5*dnorm(xx,14.5,0.75)
# lines(xx, y, col='red')

pdf2 = function(n, mu, sigma, p)
{
  # A generalization for sampling from a GMM:
  # g(x) = p_1 * N(mu_1, sigma_1) + p_2 * N(mu_2, sigma_2) + ...
  
  # assume p, mu and sigma have the same length
  n_mixt = length(p)
  
  idx = sample(1:n_mixt, size=n, replace=TRUE, prob=p)
  x = rnorm(n, mean=mu[idx], sd=sigma[idx])
  
  return (x)
}
# x = pdf2(1000, mu=c(0, 6.5,14.5), sigma=c(0.5,1.25,0.75), p=c(0.2, 0.3, 0.5))

pdf3 = function(n, mu, sigma, p)
{
  # A generalization for sampling from a GMM:
  # g(x) = p_1 * N(mu_1, sigma_1) + p_2 * N(mu_2, sigma_2) + ...
  
  # assume p, mu and sigma have the same length
  n_mixt = length(p)
  
  idx = sample(1:n_mixt, size=n, replace=TRUE, prob=p)
  x = rnorm(n)*sigma[idx] + mu[idx]
  
  return (x)
}
# x = pdf3(1000, mu=c(0, 6.5,14.5), sigma=c(0.5,1.25,0.75), p=c(0.2, 0.3, 0.5))

pdf4 = function(n, mu, sigma, p)
{
  # A generalization for sampling from a GMM:
  # g(x) = p_1 * N(mu_1, sigma_1) + p_2 * N(mu_2, sigma_2) + ...
  
  # assume p, mu and sigma have the same length
  n_mixt = length(p)
  
  idx = findInterval(runif(n), cumsum(p)) + 1
  x = rnorm(n)*sigma[idx] + mu[idx]
  
  return (x)
}
# x = pdf4(1000, mu=c(0, 6.5,14.5), sigma=c(0.5,1.25,0.75), p=c(0.2, 0.3, 0.5))