Manually generate random sample in Gaussian mixture model I want to generate (manually) a random sample in the Gaussian mixture model:
$$f_{\theta}(x) = \sum_{k = 1}^{K}\pi_k f_{\mathcal N(\mu_k, \sigma^2_k)}(x)$$
Here is my work:
theta = list(pi = c(p1,...,pK), mu = c(m1,...,mK), sigma = c(s1,...,sK))

rnormmix = function(n,theta){
  x = runif(n)

  # find mixed cdf
  pnormmix = function(x,theta){
    cdf = sapply(x, function(y) sum(theta$pi*pnorm(y,theta$mu, theta$sigma)))
    return(cdf)
  }

  # find the root of the equation fmix = F(x) - u
  fmix = function(x,u) pnormmix(x,theta) - u
  my_root = function(x) uniroot(fmix, c(-1000, 1000), tol = 0.0001, u = x)$root

  # generate sample
  return(sapply(x, my_root))
}

My questions are:


*

*Is my work correct? In the uniroot function, I impose the root must be in the interval c(-1000,1000). What happens if the root falls outside this interval?

*Are there any better methods?
 A: In principle, this seems fine, provided there is negligible mass outside of $(-1000,1000)$. It also seems problematic to extend this to a multivariate context.
A common, and more general method is to draw an assignment variable, $y \sim \operatorname{Multinomial}(\pi)$, and condition on it.
In your case this would entail:
Step 1: $y_i \stackrel{iid}{\sim} \operatorname{Multinomial}(\pi)$
Step 2: $x_i | y_i \sim \operatorname{N}(\mu_{y_i}, \sigma^2_{y_i})$
Edit: Example R code
K <- 10
pi_raw <- rexp(K)
theta <- list(pi = pi_raw/sum(pi_raw),
             mu = rnorm(K, 0, 10),
             sigma = abs(rnorm(K, 0, 3)))

rnormmix <- function(n,theta){
  y <-apply(rmultinom(n, size = 1, prob = theta$pi),2,function(col) which(col==1))
  # generate sample
  return(data.frame(x=rnorm(n, theta$mu[y], theta$sigma[y]),y=y))
}

df <- rnormmix(1000, theta)
library(ggplot2)
ggplot(data=df, aes(x=x)) + geom_histogram(aes(y=..density..)) +
  geom_segment(data=data.frame(x=theta$mu, yend=theta$pi),
               aes(x=x,xend=x,y=0,yend=yend), color="red")


A: Pick $k$ at random according to $\pi$ and simulate the corresponding normal distribution.
rnormix <- function(n, pi, mu, sigma){
  K <- length(pi)
  simulations <- numeric(n)
  for(i in 1:n){
    k <- sample.int(K, 1L, prob=pi)
    simulations[i] <- rnorm(1, mu[k], sigma[k])
  }
  return(simulations)
}

Example:
pi <- c(1/4, 3/4)
mu <- c(0,10)
sigma <- c(1,1)
sims <- rnormix(1000, pi, mu, sigma)
plot(density(sims))


For a better performance:
rnormix <- function(n, pi, mu, sigma){
  K <- length(pi)
  k <- sample.int(K, n, replace=TRUE, prob=pi)
  simulations <- rnorm(n, mu[k], sigma[k])
  return(simulations)
}

