R : function to generate a mixture distribution

Raaj picture Raaj · May 5, 2014 · Viewed 8.8k times · Source

I need to generate samples from a mixed distribution

  • 40% samples come from Gaussian(mean=2,sd=8)

  • 20% samples come from Cauchy(location=25,scale=2)

  • 40% samples come from Gaussian(mean = 10, sd=6)

To do this, i wrote the following function :

dmix <- function(x){
prob <- (0.4 * dnorm(x,mean=2,sd=8)) + (0.2 * dcauchy(x,location=25,scale=2)) + (0.4 * dnorm(x,mean=10,sd=6))
return (prob)
}

And then tested with:

foo = seq(-5,5,by = 0.01)
vector = NULL
for (i in 1:1000){
vector[i] <- dmix(foo[i])
}
hist(vector)

I'm getting a histogram like this (which I know is wrong) - histogram of the supposed distribution

What am I doing wrong? Can anyone give some pointers please?

Answer

Josh O&#39;Brien picture Josh O'Brien · May 5, 2014

There are of course other ways to do this, but the distr package makes it pretty darned simple. (See also this answer for another example and some more details about distr and friends).

library(distr)

## Construct the distribution object.
myMix <- UnivarMixingDistribution(Norm(mean=2, sd=8), 
                                  Cauchy(location=25, scale=2),
                                  Norm(mean=10, sd=6),
                                  mixCoeff=c(0.4, 0.2, 0.4))
## ... and then a function for sampling random variates from it
rmyMix <- r(myMix)

## Sample a million random variates, and plot (part of) their histogram
x <- rmyMix(1e6)
hist(x[x>-100 & x<100], breaks=100, col="grey", main="")

enter image description here

And if you'd just like a direct look at your mixture distribution's pdf, do:

plot(myMix, to.draw.arg="d") 

enter image description here