[R] nlminb supplying NaN parameters to objective function

[Jean Marchal]

Dear list,

I am doing some maximum likelihood estimation using nlminb() with
box-constraints to ensure that all parameters are positive. However,
nlminb() is behaving strangely and seems to supply NaN as parameters
to my objective function (confirmed using browser()) and output the
following:

$par
 [1] NaN NaN NaN   0 NaN   0 NaN NaN NaN NaN NaN NaN NaN

$objective
[1] 0

$convergence
[1] 1

$iterations
[1] 19

$evaluations
function gradient
      87      542

$message
[1] "gr cannot be computed at initial par (65)"


When I use trace = TRUE, I can see the following:

  0:     32495.488: 0.0917404 0.703453  1.89661 1.11022e-16
1.11022e-16 0.107479 1.11022e-16 1.11022e-16 1.11022e-16 0.472377
0.894128  1.86743 1.11022e-16
  1:     4035.3900: 0.0917404 0.703453  1.89661 1.11022e-16
1.11022e-16 0.107479 1.11022e-16 1.11022e-16 1.11022e-16 0.472377
0.894128  1.86743 0.250000
  2:     3955.8801: 0.0948452 0.704168  1.89651 0.000135456 0.0310485
0.107991 0.00138902 0.000427631 1.11022e-16 0.472331 0.894128  1.86743
0.250000
  3:     3951.4141: 0.0948926 0.703906  1.89640 2.99167e-05 0.0315288
0.109692 0.00242572 0.00272185 7.96814e-05 0.472780 0.894130  1.86744
0.249998
....
 17:     3937.3923: 0.0947470 0.703030  1.89605 1.11022e-16 0.0300763
0.115081 0.00562496 0.00989997 0.000323268 0.474247 0.894142  1.86745
0.249737
 18:     3937.3923: 0.0947470 0.703030  1.89605 1.11022e-16 0.0300763
0.115081 0.00562496 0.00989997 0.000323268 0.474247 0.894142  1.86745
0.249737
 19:    -0.0000000:     -nan     -nan     -nan 1.11022e-16     -nan
 -nan     -nan     -nan     -nan     -nan     -nan     -nan      nan


my objective function looks like:

nLL <- function(params){

  mu <- drop(model.matrix(modelTermsObj) %*% params)

  if(any(mu < 0) || anyNA(mu) || any(is.infinite(mu))){
    return(.Machine$double.xmax)
  } else {
    return(-sum(dnbinom(x=args$data[,response], mu = mu, size =
params[length(params)], log = TRUE)))
  }
}

I tried different starting values, different bounds but without
success so far. Is this a bug?

PS after trying to make a reproducible example that I gracefully
failed to do... I change my objective function so instead of using
model.matrix(), I did the maths (e.g. Y ~ A + B * C). Thus, mu is now
a bunch of NaN, and my objective function return .Machine$double.xmax
which is fine. Then nlminb stops and returns (like if nothing
happened):

$par
 [1] 1.11022e-16 1.11022e-16 2.69205e-04 1.11022e-16 1.68161e-03
1.06027e-03 1.16969e-05 1.11022e-16 8.51669e+01 7.31162e+01
5.04748e+00 5.28373e+00 1.23992e-01

$objective
[1] 3823.567

$convergence
[1] 0

$iterations
[1] 1

$evaluations
function gradient
       2       13

$message
[1] "X-convergence (3)"

I can provide the data and model if necessary but cannot make them
publicly available (yet).

Thank you,

Jean