Re: [R] Problem with mle

From: Ben Bolker <bolker_at_ufl.edu>
Date: Sat 03 Jun 2006 - 03:42:21 EST

Rainer M KRug <RMK <at> krugs.de> writes:

>
> R 2.3.0
> Linux, SuSE 10.0
>
> Hi
>
> I have two problems with mle - probably I am using it the wrong way so
> please let me know.
>
> I want to fit different distributions to an observed count of seeds and
> in the next step use AIC or BIC to identify the best distribution.

library(stats4)
SumSeeds <- c(762, 406, 288, 260, 153, 116, 57, 33, 40, 44, 36,

              24, 35, 23, 36, 25, 35, 30) X <- c(1.74924, 3.49848, 5.24772, 6.99696, 8.74620, 17.49240, 26.23860,

       34.98480, 43.73100, 52.47720, 61.22340, 69.96960, 71.71880, 73.46810,
       75.21730, 76.96650, 78.71580, 80.46500 )

plot(X,SumSeeds)
curve(7000*dexp(x,0.1),add=TRUE,from=0)

## exponential function (a*exp(-rate*X)) makes more sense
##  than using dexp, which is normalized to 1 (thus makes
## it harder to interpret a as an intercept -- I don't
## know what X is?  (You could easily reverse this decision
## though.)

negexplik <- function(a=7000,rate=0.1) {   -sum(dpois(SumSeeds,lambda=a*dexp(X,rate),log=TRUE)) }

## start from eyeball-estimated parameters m1 = mle(minuslogl=negexplik,start=list(a=7000,rate=0.1)) ## add estimated curve
with(as.list(coef(m1)),curve(a*dexp(x,rate),col="red",add=TRUE,from=0))

## conclusion so far: the exponential fit is really bad, because ## the numbers don't fall to zero as fast as expected

## changed parameterization of Weibull slightly to agree with ## the parameterization of the exponential above weiblik <- function(a=7000,shape=1,rate=0.1) {   -sum(dpois(SumSeeds,lambda=a*dweibull(X,scale=1/rate,shape=shape),

        log=TRUE))
}

weiblik(a=7222,shape=1,rate=0.05) ## check ## start from exponential a/rate estimates, plus shape=1 which reduces ## to exponential
m2 = mle(minuslogl=weiblik,start=list(a=7222,rate=0.05,shape=1)) with(as.list(coef(m2)),curve(a*dweibull(x,scale=1/rate,shape=shape),

     col="blue",add=TRUE,from=0))

as.numeric(pchisq(-2*(logLik(m1)-logLik(m2)),df=1,

      lower.tail=FALSE,log.p=TRUE))
## insanely significant

  I still had some problems profiling. Various things that may help: (1) set parscale= option
(2) do fits on the log-parameters (which all have to be positive) OR   (3) use L-BFGS-B and set lower=0
(4) summary(m1) will give you approximate (Fisher information) confidence limits even if profiling doesn't work.

  cheers
    Ben Bolker



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Jun 03 04:56:47 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 05 Jun 2006 - 20:10:32 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.