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

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