20
$\begingroup$

How can I sample from a mixture distribution, and in particular a mixture of Normal distributions in R? For example, if I wanted to sample from:

$$ 0.3\!\times\mathcal{N}(0,1)\; + \;0.5\!\times\mathcal{N}(10,1)\; + \;0.2\!\times\mathcal{N}(3,.1) $$

how could I do that?

$\endgroup$
  • 3
    $\begingroup$ I really don't like this way of denoting a mixture. I know it's conventionally done like this but I find it misleading.The notation suggests that to sample, you need to sample all three normals and weigh the results by those coefficients which would obviously not be correct. Anyone know a better notation? $\endgroup$ – StijnDeVuyst May 28 '16 at 19:14
  • $\begingroup$ I never got that impression. I think of the distributions (in this case the three normal distributions) as functions and then the result is another function. $\endgroup$ – roundsquare Oct 12 at 20:13
  • $\begingroup$ @StijnDeVuyst you might want to visit this question originated from your comment : stats.stackexchange.com/questions/431171/… $\endgroup$ – ankii Oct 12 at 20:25
  • $\begingroup$ @ankii: thanks for pointing that out! $\endgroup$ – StijnDeVuyst Oct 14 at 15:51
32
$\begingroup$

It's good practice to avoid for loops in R for performance reasons. An alternative solution which exploits the fact rnorm is vectorized:

N <- 100000

components <- sample(1:3,prob=c(0.3,0.5,0.2),size=N,replace=TRUE)
mus <- c(0,10,3)
sds <- sqrt(c(1,1,0.1))

samples <- rnorm(n=N,mean=mus[components],sd=sds[components])
$\endgroup$
  • 3
    $\begingroup$ Alternatively, you can use the properties of the normal distribution to replace the last line by samples <- rnorm(N)*sds[components]+mus[components]. I find it easier to read :) $\endgroup$ – Elvis Sep 24 '13 at 8:59
  • $\begingroup$ Very elegant (cc @Elvis) ! $\endgroup$ – Itamar Oct 2 '13 at 7:13
18
$\begingroup$

In general, one of the easiest ways to sample from a mixture distribution is the following:

Algorithm Steps

1) Generate a random variable $U\sim\text{Uniform}(0,1)$

2) If $U\in\left[\sum_{i=1}^kp_{k},\sum_{i=1}^{k+1}p_{k+1}\right)$ interval, where $p_{k}$ correspond to the the probability of the $k^{th}$ component of the mixture model, then generate from thedistribution of the $k^{th}$ component

3) Repeat steps 1) and 2) until you have the desired amount of samples from the mixture distribution

Now using the general algorithm given above, you could sample from your example mixture of normals by using the following R code:

#The number of samples from the mixture distribution
N = 100000                 

#Sample N random uniforms U
U =runif(N)

#Variable to store the samples from the mixture distribution                                             
rand.samples = rep(NA,N)

#Sampling from the mixture
for(i in 1:N){
    if(U[i]<.3){
        rand.samples[i] = rnorm(1,0,1)
    }else if(U[i]<.8){
        rand.samples[i] = rnorm(1,10,1)
    }else{
        rand.samples[i] = rnorm(1,3,.1)
    }
}

#Density plot of the random samples
plot(density(rand.samples),main="Density Estimate of the Mixture Model")

#Plotting the true density as a sanity check
x = seq(-20,20,.1)
truth = .3*dnorm(x,0,1) + .5*dnorm(x,10,1) + .2*dnorm(x,3,.1)
plot(density(rand.samples),main="Density Estimate of the Mixture Model",ylim=c(0,.2),lwd=2)
lines(x,truth,col="red",lwd=2)

legend("topleft",c("True Density","Estimated Density"),col=c("red","black"),lwd=2)

Which generates:

enter image description here

and as a sanity check:

enter image description here

$\endgroup$
  • $\begingroup$ Hi! Thanks so much! This answer helped me greatly. I am using this in a research project. I wish to quote a reference for the above. Can you please suggest a research article citation. $\endgroup$ – Abhishek Bhatia Jun 29 '15 at 11:24
7
$\begingroup$

Conceptually, you are just picking one distribution (from $k$ possibilities) with some probability, and then generating pseudo-random variates from that distribution. In R, this would be (e.g.):

set.seed(8)               # this makes the example reproducible
N     = 1000              # this is how many data you want
probs = c(.3,.8)          # these are *cumulative* probabilities; since they 
                          #   necessarily sum to 1, the last would be redundant
dists = runif(N)          # here I'm generating random variates from a uniform
                          #   to select the relevant distribution

# this is where the actual data are generated, it's just some if->then
#   statements, followed by the normal distributions you were interested in
data = vector(length=N)
for(i in 1:N){
  if(dists[i]<probs[1]){
    data[i] = rnorm(1, mean=0, sd=1)
  } else if(dists[i]<probs[2]){
    data[i] = rnorm(1, mean=10, sd=1)
  } else {
    data[i] = rnorm(1, mean=3, sd=.1)
  }
}

# here are a couple of ways of looking at the results
summary(data)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# -3.2820  0.8443  3.1910  5.5350 10.0700 13.1600 

plot(density(data))

enter image description here

$\endgroup$
  • $\begingroup$ Nice answer, you beat me to posting :P $\endgroup$ – user25658 Sep 24 '13 at 1:23
  • 1
    $\begingroup$ Thanks for the tip, @BabakP. I'm not sure what it was. It was something in the ifelse() statement, but I'll have to figure it out later. I replaced that code w/ a loop. $\endgroup$ – gung - Reinstate Monica Sep 24 '13 at 2:32
  • 6
    $\begingroup$ (cc @BabakP) These are both good answers and are obviously correct (+1s). Just an R programming trick: you can also use the findInterval() and cumsum() commands to simplify the code and, more importantly, make it easier to generalize to a different number of dimensions. For example, for an input vector of means $\mu$ (mu) and variances $\sigma^2$(s), and mixture probabilities (p), a simple function to generate n samples from this mixture would be mix <- function(n,mu,s,p) { ii <- findInterval(runif(n),cumsum(p))+1; x <- rnorm(n,mean=mu[ii],sd=sqrt(s[ii])); return(x); } $\endgroup$ – Macro Sep 24 '13 at 3:19
  • 1
    $\begingroup$ @Macro, very true and very nice code! I have not seen the findInterval() command before, however, I like to write code on here as simplistically as I can because I want it to be a tool for understanding rather than efficiency. $\endgroup$ – user25658 Sep 24 '13 at 3:22
  • 1
    $\begingroup$ I said these were good answers. My purpose was not to criticize you but to offer an approach that easily generalizes to more than three dimensions by only changing a single argument, not any code. It's not clear to me why what you've written is any more transparent than what I've written but I surely don't want to argue about that. Cheers. $\endgroup$ – Macro Sep 24 '13 at 3:54
0
$\begingroup$

Already given perfect answers, so for those who want to achieve this in Python, here is my solution:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

mu = [0, 10, 3]
sigma = [1, 1, 1]
p_i = [0.3, 0.5, 0.2]
n = 10000

x = []
for i in range(n):
    z_i = np.argmax(np.random.multinomial(1, p_i))
    x_i = np.random.normal(mu[z_i], sigma[z_i])
    x.append(x_i)

def univariate_normal(x, mean, variance):
    """pdf of the univariate normal distribution."""
    return ((1. / np.sqrt(2 * np.pi * variance)) * 
            np.exp(-(x - mean)**2 / (2 * variance)))

a = np.arange(-7, 18, 0.01)
y = p_i[0] * univariate_normal(a, mean=mu[0], variance=sigma[0]**2) + p_i[1] * univariate_normal(a, mean=mu[1], variance=sigma[0]**2)+ p_i[2] * univariate_normal(a, mean=mu[2], variance=sigma[0]**2)

fig, ax = plt.subplots(figsize=(8, 4))

ax.hist(x, bins=100, density=True)
ax.plot(a, y)

enter image description here

$\endgroup$

Your Answer

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