The problem is that fitting mixtures models is hard -- (non)identifiability is a serious problem: very large sample sizes may be necessary to clearly distinguish the modes. As V&R say in MASS, 4th edition, p. 320: " ... fitting normal mixtures is a difficult problem, and the results obtained are often heavily dependent on the initial configuration ([i.e. starting values. BG] supplied"
The EM algorithm is quite commonly used: you might have a look at em() and related functions in the mclust package if Ben's straightforward approach fails to do the job for you. No guarantee em will work either, of course. -- Bert Bert Gunter Genentech Nonclinical Biostatistics 650-467-7374 -----Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Ben Bolker Sent: Wednesday, April 08, 2009 12:47 PM To: r-help@r-project.org Subject: Re: [R] MLE for bimodal distribution _nico_ wrote: > > Hello everyone, > > I'm trying to use mle from package stats4 to fit a bi/multi-modal > distribution to some data, but I have some problems with it. > Here's what I'm doing (for a bimodal distribution): > > # Build some fake binormally distributed data, the procedure fails also > with real data, so the problem isn't here > data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3)) > # Just to check it's bimodal > plot(density(data)) > f = function(m, s, m2, s2, w) > { > -log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2, > sd=s2)) ) > } > > res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6)) > > This gives an error: > Error in optim(start, f, method = method, hessian = TRUE, ...) : > non-finite finite-difference value [2] > In addition: There were 50 or more warnings (use warnings() to see the > first 50) > And the warnings are about dnorm producing NaNs > > So, my questions are: > 1) What does "non-finite finite-difference value" mean? My guess would be > that an Inf was passed somewhere where a finite number was expected...? > I'm not sure how optim works though... > 2) Is the log likelihood function I wrote correct? Can I use the same type > of function for 3 modes? > 3) What should I do to avoid function failure? I tried by changing the > parameters, but it did not work. > 4) Can I put constraints to the parameters? In particular, I would like w > to be between 0 and 1. > 5) Is there a way to get the log-likelihood value, so that I can compare > different extimations? > 6) Do you have any (possibly a bit "pratically oriented") readings about > MLE to suggest? > > Thanks in advance > Nico > Here's some tweaked code that works. Read about the "L-BFGS-B" method in ?optim to constrain parameters, although you will have some difficulty making this work for more than two components. I think there's also a mixture model problem in Venables & Ripley (MASS). There are almost certainly more efficient approaches for mixture model fitting, although this "brute force" approach should be fine if you don't need to do anything too complicated. (RSiteSearch("mixture model")) set.seed(1001) data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3)) # Just to check it's bimodal plot(density(data)) f = function(m, s, m2, s2, w) { -sum(log(w*dnorm(data, mean=m, sd=s)+ (1-w)*dnorm(data, mean=m2, sd=s2))) } library(stats4) start0 <- list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6) do.call("f",start0) ## OK res = mle(f, start=start0) with(as.list(coef(res)), curve(w*dnorm(x,m,s)+(1-w)*dnorm(x,m2,s2),add=TRUE,col=2)) -- View this message in context: http://www.nabble.com/MLE-for-bimodal-distribution-tp22954970p22957613.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.