I would like to fit a mixture model to Monte Carlo generated data with probability densities which typically look like those in the attached image. typical densities

It would seem from visual inspection that a normal mixture model might be applicable but on checking the CRAN task view I really do not know which package might be appropriate for my needs.

Basically what I would like to do is supply a vector of the data and then have the package function return the mean, variance and proportional weights for each component in the mixture model, and also perhaps identify how many components there are in the model.


Try mixdist

Here's an example:


#Build data vector "x" as a mixture of data from 3 Normal Distributions  
x1 <- rnorm(1000, mean=0, sd=2.0)  
x2 <- rnorm(500, mean=9, sd=1.5)  
x3 <- rnorm(300, mean=13, sd=1.0)  
x <- c(x1, x2, x3)  

#Plot a histogram (you'll play around with the value for "breaks" as    
#you zero-in on the fit).   Then build a data frame that has the  
#bucket midpoints and counts.  
breaks <- 30  
his <- hist(x, breaks=breaks)  
df <- data.frame(mid=his$mids, cou=his$counts)  

#The above Histogram shows 3 peaks that might be represented by 3 Normal  
#Distributions.  Guess at the 3 Means in Ascending Order, with a guess for  
#the associated 3 Sigmas and fit the distribution.  
guemea <- c(3, 11, 14)  
guesig <- c(1, 1, 1)  
guedis <- "norm"  
(fitpro <- mix(as.mixdata(df), mixparam(mu=guemea, sigma=guesig), dist=guedis))  

#Plot the results  
plot(fitpro, main="Fit a Probability Distribution")  
legend("topright", lty=1, lwd=c(1, 1, 2), c("Original Distribution to be Fit", "Individual Fitted Distributions", "Fitted Distributions Combined"), col=c("blue", "red", rgb(0.2, 0.7, 0.2)), bg="white")  


      pi     mu  sigma  
1 0.5533 -0.565 1.9671  
2 0.2907  8.570 1.6169  
3 0.1561 12.725 0.9987  

[1] "norm"  

   conpi    conmu consigma   
  "NONE"   "NONE"   "NONE"   

enter image description here


Package Mclust is nice. The mclust function fits a mixture of normals distribution to data. You can automatically choose the number of components based on BIC (mclustmodel) or specify the number of components. There is also no need to convert your data into a data frame.

Also, package Mixtools and the function normalmixEM fits a mixture of normals.

Update: I recently discovered the mixAK package and the NMixMCMC function and it is terrific. It has many options, including RJMCMC for component selection, right left censoring, etc...


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.