Installation of multimodal distributions in R; generating new values ​​from the established distribution

Question

I work 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 is as follows:

I know that the values ​​are taken from two modes - low and high - and if the main process is normal, I used the mixtools package to match the 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:

(solid lines are curves, and a 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 

Mix Parameters:

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

• To create a new set of values ​​that approaches the base distribution, is my approach correct or is there a better workflow?
• 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 little dubious to be suitable mixes, but ignore it. You can take a sample from the installed 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]) 

Answer 2

Your approach is right.

For each sample from your mixed distribution, you just need to choose which of the two component Gaussian distributions the sample should take, and then extract the sample from this distribution.

You can choose between two distributions using the proportions of the mixture you found: simulate a random number between 0 and 1 and a sample from the first distribution, if this random number is less than the first proportion, otherwise the sample from the second distribution.

Finally, a sample from the corresponding Gaussian distribution using the rnorm function.

 dat.demand2050.mixmdl.prop=c(0.2783939,0.7216061) dat.demand2050.mixmdl.means=c(56.21150,73.08389) dat.demand2050.mixmdl.dev=c(3.098292,6.413906) sampleMixture=function(prop,means,dev){ # Generate a uniformly distributed random number between 0 and 1 # in order to choose between the two component distributions distTest=runif(1) if(distTest<prop[1]){ # Then sample from the first component of the mixture sample=rnorm(1,mean=means[1],sd=dev[1]) }else{ # Sample from the second component of the mixture sample=rnorm(1,mean=means[2],sd=dev[2]) } return(sample) } # Generate a single sample sampleMixture(dat.demand2050.mixmdl.prop,dat.demand2050.mixmdl.means,dat.demand2050.mixmdl.dev) # Generate 100 samples and plot resulting distribution samples=replicate(100,sampleMixture(dat.demand2050.mixmdl.prop,dat.demand2050.mixmdl.means,dat.demand2050.mixmdl.dev)) plot(density(samples))