[R] R 3.1.2 mle2() function on Windows 7 Error and multiple optimal solutions

[Ben Bolker]

Mosab Alqurashi <malqurashi <at> KSU.EDU.SA> writes:

> 
> ?Hello,
 
> I want to get the Maximum Likelihood Estimate (MLE) for an Extended
> weibull with three parameters. I have two problems. First, some
> times if I change the starting point for one of the parameters I got
> an error message. Second, sometime when I change the starting point
> for one of the parameters the solution changes. I don't know what
> I'm doing wrong or if there is a better way to get the MLE for a
> distribution.

  I haven't had a chance to look closely at this, but some possibilities
are:
 
  * there are actually multiple optima for this function (test: if
you start at or near one of the optima, do you stay/return to the
same optimum?)

 * some of your starting values might be very far away from the
true values, in which case the likelihood surface could be very flat

 * some of your parameters could be very strongly correlated,
in which case the likelihood surface is very flat in some directions

 * if the parameters are on very different scales you might try
adjusting the "parscale" option (see ?optim)

 * it might be worth experimenting with other optimizers (although
L-BFGS-B is the only optimizer in base R that allows constraints)

 * you might try fitting your parameters on the log scale, to 
avoid needing the constraints

 * you could try AD Model Builder 

 * you could try visualizing the likelihood surface to see what's
going on
 
> Here is my Example and the Error message.
> 
> Problem 1:
> 
> y=c(5,11,21,31,46,75,98,122,145,165,195,224,245,293,321,330,350,420) 
> #My data set
> 
> Nweibull <-
> function(x,lambda,alpha,beta){
    lambda*beta*(x/alpha)^(beta-1)*exp(1)^((x/alpha)^beta+
      lambda*alpha*(1-exp(1)^((x/alpha)^beta)))}
>  # The PDF
> 
> LL <- function(lambda,alpha,beta) {
>   R = Nweibull(y,lambda,alpha,beta)
> 
>   -sum(log(R))
> } # defining the Likelihood Function
> 
> g=mle2(LL, start = list(lambda=.01,alpha=325,beta=.8),
> data = list(y),method =
> "L-BFGS-B",lower=c(.00001,.00001,.00001),upper = c(Inf, Inf,Inf))
> summary(g)
> 
> #### Out Put###
> 
> Maximum likelihood estimation
> Call:
> mle2(minuslogl = LL, start = list(lambda = 0.01, alpha = 325,
>     beta = 0.8), method = "L-BFGS-B", data = list(y), lower = c(1e-05,
>     1e-05, 1e-05), upper = c(Inf, Inf, Inf))
> Coefficients:
>          Estimate Std. Error    z value     Pr(z)
> lambda 3.6926e-03 8.9662e-04 4.1183e+00 3.817e-05 ***
> alpha  3.2500e+02 1.1953e-04 2.7191e+06 < 2.2e-16 ***
> beta   9.3856e-01 2.0186e-01 4.6496e+00 3.326e-06 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> -2 log L: 218.7464
> 
> Now If I change the starting point for theta to be 0.7 
> I got the following error message:
> g=mle2(LL, start = list(lambda=.01,alpha=340,beta=.7),data = list(y),method =
> "L-BFGS-B",lower=c(.00001,.00001,.00001),upper = c(Inf, Inf,Inf))
> 
> Error in optim(par = c(0.01, 325, 0.7), fn = function (p)  :
>   L-BFGS-B needs finite values of 'fn'
> 
> Problem 2:
> If I change the starting point for alpha to be 
> 340 the optimal solution will change too:
> g=mle2(LL, start = list(lambda=.01,alpha=340,beta=.8),data = list(y),method =
> "L-BFGS-B",lower=c(.00001,.00001,.00001),upper = c(Inf, Inf,Inf))
> summary(g)
> 
> Maximum likelihood estimation
> 
> Call:
> mle2(minuslogl = LL, start = list(lambda = 0.01, alpha = 340,
>     beta = 0.8), method = "L-BFGS-B", data = list(y), lower = c(1e-05,
>     1e-05, 1e-05), upper = c(Inf, Inf, Inf))
> 
> Coefficients:
>          Estimate Std. Error    z value     Pr(z)
> lambda 3.7938e-03 9.3580e-04 4.0541e+00 5.034e-05 ***
> alpha  3.4000e+02 1.2039e-04 2.8243e+06 < 2.2e-16 ***
> beta   9.6164e-01 2.0766e-01 4.6309e+00 3.642e-06 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> -2 log L: 218.7979
> 
> Thanks,
> Mosab
>