[R] MLE for bimodal distribution

[Ben Bolker]

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