Fitting multimodal distributions in R; generating new values from fitted distribution

Question

I am working with small sample size data:

>dput(dat.demand2050.unique)  
c(79, 56, 69, 61, 53, 73, 72, 86, 75, 68, 74.2, 80, 65.6, 60, 54)    

for which the density distribution looks like this:

I know that the values are from two regimes - low and high - and assuming that the underlying process is normal, I used the mixtools package to fit a bimodal distribution:

set.seed(99)  
dat.demand2050.mixmdl <- normalmixEM(dat.demand2050.unique, lambda=c(0.3,0.7), mu=c(60,70), k=2)

which gives me the following result:

(the solid lines are fitted curves and the dashed line is the original density).

# get the parameters of the mixture
dat.demand2050.mixmdl.prop <- dat.demand2050.mixmdl$lambda    #mix proportions
dat.demand2050.mixmdl.means <- dat.demand2050.mixmdl$mu    #modal means
dat.demand2050.mixmdl.dev <- dat.demand2050.mixmdl$sigma   #modal std dev  

The mixture parameters are:

>dat.demand2050.mixmdl.prop  #mix proportions  
[1] 0.2783939 0.7216061  
>dat.demand2050.mixmdl.means  #modal means  
[1] 56.21150 73.08389  
>dat.demand2050.mixmdl.dev  #modal std dev  
[1] 3.098292 6.413906 

I have the following questions:

1. To generate a new set of values that approximates the underlying distribution, is my approach correct or is there a better workflow?
2. If my approach is correct, how can I use this result to generate a set of random values from this mixed distribution?

Answer 1

Your sample size is a bit dubious to be fitting mixtures, but never mind that. You can sample from the fitted mixture as follows:

probs <- dat.demand2050.mixmdl$lambda
m <- dat.demand2050.mixmdl$mu
s <- at.demand2050.mixmdl$sigma

N <- 1e5
grp <- sample(length(probs), N, replace=TRUE, prob=probs)
x <- rnorm(N, m[grp], s[grp])