How do I find best starting values for nlimb optimization?

Question

I am trying to run the following code with the attached dataset. How do I solve the error of hessian matrix inversion?

library(stats4)
library(bbmle)
library(stats)
library(numDeriv)
library('bbmle')
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3, 1.7, 2.3, 1.6, 2.0)
hist(x)
fEHLKUMW<-function(a,b,alpha,vartheta)
{
  -sum(log(  (2*a*b*alpha*vartheta*(x^(vartheta-1))*(exp(-x^(vartheta)))*((1- (exp(-x^(vartheta))))^(a-1))*((1-((1-((1-(exp(-x^(vartheta))))^a))^b))^(alpha-1)))/((1-((1-(exp(-x^(vartheta))))^a))^(b*(alpha+1)))
))
}
EHLKUMW.result<-mle2(fEHLKUMW,hessian = NULL,start=list(a=0.01,b=0.01,alpha=.3,vartheta=0.01),optimizer="nlminb",lower=0)
summary(EHLKUMW.result)

I get the error as;

**
Warning messages:
1: In nlminb(start = start, objective = objectivefunction, hessian = NULL,  :
  NA/NaN function evaluation
2: In mle2(fEHLKUML, hessian = NULL, start = list(a = 1, b = 0.4, c = 0.5,  :
  couldn't invert Hessian
**

Answer 1

This is a very open question I think, so I'll present some tools and an approach to this, and maybe others can comment, etc.

First, the main part:

library(bbmle)
library(stats)
library(numDeriv)
library(bbmle)
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3, 1.7, 2.3, 1.6, 2.0)
hist(x)
fEHLKUMW <- function(a,b,alpha,vartheta) {
  -sum(log(  (2*a*b*alpha*vartheta*(x^(vartheta-1))*(exp(-x^(vartheta)))*((1- (exp(-x^(vartheta))))^(a-1))*((1-((1-((1-(exp(-x^(vartheta))))^a))^b))^(alpha-1)))/((1-((1-(exp(-x^(vartheta))))^a))^(b*(alpha+1)))
  ))
}

Now, we can of course run it like you've done:

EHLKUMW.result <- mle2(
  fEHLKUMW,
  hessian = NULL,
  start = list(
    a = 0.01,
    b = 0.01,
    alpha = .3,
    vartheta = 0.01
  ),
  optimizer = "nlminb",
  lower = 0
)

But we can also run it with a distribution on each of these parameters, to get a new input all the time:

EHLKUMW.result <- mle2(
  fEHLKUMW,
  hessian = NULL,
  start = 
    list(
      # a = 0.01,
      # a = rt(1, 10, ncp = 0.01),
      a = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
      # b = 0.01,
      b = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
      # alpha = .3,
      alpha = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.3, df = 10),
      # vartheta = 0.01
      vartheta = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10)
    ),    
  optimizer = "nlminb",
  lower = 0
)

I've chosen to use trundist to get a t-distribution, centralised around the values you provided, and lower is 0 through the a = -argument. If you know what the upper limit of these parameters, this can be done with b = -argument.

The output that I think is most relevant are the attained logLik and the coef.

library(tidyverse)

exec(fEHLKUMW, !!!list(
  a = 0.01,
  b = 0.01,
  alpha = .3,
  vartheta = 0.01
))
replicate(
  250,
  exec(fEHLKUMW, !!!list(
    # a = 0.01,
    # a = rt(1, 10, ncp = 0.01),
    a = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
    # b = 0.01,
    b = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
    # alpha = .3,
    alpha = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.3, df = 10),
    # vartheta = 0.01
    vartheta = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10)
  )))

I used this to fine-tune the distributions the parameters

The following is a multiple runs and their output compared to logLik.

tibble(n = seq_len(100),
       output = map(n, ~mle2(
         fEHLKUMW,
         hessian = NULL,
         start = 
           list(
             # a = 0.01,
             # a = rt(1, 10, ncp = 0.01),
             a = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
             # b = 0.01,
             b = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
             # alpha = .3,
             alpha = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.3, df = 10),
             # vartheta = 0.01
             vartheta = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10)
           ),    
         optimizer = "nlminb",
         lower = 0
       ))) ->
  outputs_df

This code gives a nice print

outputs_df %>% 
  mutate(coef = output %>% map(coef),
         logLik = output %>% map_dbl(logLik)) %>% 
  unnest_wider(coef) %>% 
  arrange(logLik) %>% 
  print(n=Inf)