How to insert a numeric gradient inside constrOptim

Question

I was doing maximum likelihood estimation using optim() and it was quite easy. It's a generalized logistic distribution with 4 parameters and a couple of restrictions, all listed in the likelihood function:

genlogis.loglikelihood <- function(param = c(sqrt(2/pi),0.5, 2, 0), x){

  if(length(param) < 3 | length(param) > 4 ){
    stop('Incorrect number of parameters: param = c(a,b,p,location)')
  }

  if(length(param) == 3){
    #warning('Location parameter is set to 0')
    location = 0
  }

  if(length(param) == 4){
    location = param[4]
  }

  a = param[1]
  b = param[2]
  p = param[3]

  if(!missing(a)){
    if(a < 0){
      stop('The argument "a" must be positive.')
    }
  }
  if(!missing(b)){
    if(b < 0){

      stop('The argument "b" must be positive.')
    }
  }
  if(!missing(p)){
    if(p < 0){
      stop('The argument "p" must be positive.')
    }
  }

  if(p == 0 && b > 0 && a > 0){
    stop('If "p" equals to 0, "b" or "a" must be 
         0 otherwise there is identifiability problem.')
  }  
  if(b == 0 && a == 0){
    stop('The distribution is not defined for "a" 
         and "b" equal to 0 simultaneously.')
  }

  z <- sum(log((a+b*(1+p)*abs((x-location))^p ) * exp(-((x-location)*(a+b*abs((x-location))^p))))) -
             sum(2*log(exp(-((x-location)*(a+b*abs((x-location))^p))) + 1))
  if(!is.finite(z)){
    z <- 1e+20
  }

  return(-z)
}

I made it's likelihood function and worked flawessly this way:

    opt <-  function(parameters, data){
            optim(par = parameters, fn = genlogis.loglikelihood, x=data,
                  lower = c(0.00001,0.00001,0.00001, -Inf), 
                  upper = c(Inf,Inf,Inf,Inf), method = 'L-BFGS-B') 
        }
opt(c(0.3, 1.01, 2.11, 3.5), faithful$eruptions)

Since this function does the gradient numerically I had not much problem.

Then I wanted to change to constrOptim() because the boundaries are actually 0 and not a small number on the first 3 parameters. But, the problem I face is that the argument grad has to be specified and I can't derive that function to give a gradient function, so I have to do it numerically as in optim(), it works if I put grad = NULL but I don't want Nelder-Mead method but BFGS.

I've tried this way but not of much sucess:

   opt2 <- function(initial, data){   
  ui <- rbind(c(1, 0, 0, 0), c(0,1,0,0), c(0,0,1,0))
  ci <- c(0,0,0)    
            constrOptim(theta = initial, f = genlogis.loglikelihood(param, x), 
                        grad = numDeriv::grad(func = function(x, param) genlogis.loglikelihood(param, x), param = theta, x = data)
                        , x = data, ui = ui, ci = ci)
        }

Answer 1

Your notation is a bit complicated, maybe that confused you.

opt2 <- function(parameters, data){
  fn = function(p) genlogis.loglikelihood(p, x = data)
  gr = function(p) numDeriv::grad(fn, p)
  ui <- rbind(c(1, 0, 0, 0), c(0,1,0,0), c(0,0,1,0))
  ci <- c(0,0,0)    
  constrOptim(theta = parameters, f = fn, grad = gr,
              ui = ui, ci = ci, method="BFGS")
}
opt2(c(0.3, 1.01, 2.11, 3.5), faithful$eruptions)