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?


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


| cite | improve this answer | |
  • $\begingroup$ Thanks! In case the variable is fixed (like mine), do you have any suggestion? $\endgroup$ – SiXUlm Apr 6 '17 at 17:15
  • $\begingroup$ @SiXUlm $\theta$ is fixed in the above example. Does this answer your question? $\endgroup$ – HStamper Apr 7 '17 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])
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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