Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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.

share|improve this question
add comment

2 Answers

up vote 7 down vote accepted

Try mixdist

Here's an example:

library(mixdist)  

#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)  
head(df)  

#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")  
grid()  
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")  


===========================  


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

Distribution:  
[1] "norm"  

Constraints:  
   conpi    conmu consigma   
  "NONE"   "NONE"   "NONE"   

enter image description here

share|improve this answer
add comment

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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