Raaj - 1 year ago 84

R Question

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

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

Answer Source

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

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

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