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)))

  # 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:

  1. 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?

  2. Are there any better methods?


2 Answers 2


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)
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")


  • $\begingroup$ Thanks! In case the variable is fixed (like mine), do you have any suggestion? $\endgroup$
    – SiXUlm
    Commented Apr 6, 2017 at 17:15
  • $\begingroup$ @SiXUlm $\theta$ is fixed in the above example. Does this answer your question? $\endgroup$
    – HStamper
    Commented Apr 7, 2017 at 20:58

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])


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

gaussian mixture

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])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.