[R] nlminb supplying NaN parameters to objective function

[Jean Marchal]

A follow-up to my yesterday's email.

I was able to make a reproducible example. All you will have to do is
load the .RData file that you can download here:
https://drive.google.com/file/d/0B0DKwRjF11x4dG1uRWhwb1pfQ2s/view?usp=sharing

and run this line of code:

nlminb(start=sv, objective = nLL, lower = 0, upper = Inf,
control=list(trace=TRUE))

which output the following:

  0:     12523.401: 0.0328502 0.0744493 0.00205298 0.0248628 0.0881807
0.0148887 0.0244485 0.0385922 0.0714495 0.0161784 0.0617551 0.0244901
0.0784038
  1:     12421.888: 0.0282245 0.0697934  0.00000 0.0199076 0.0833634
0.0101135 0.0189494 0.0336236 0.0712130 0.0160687 0.0616015 0.0244689
0.0660129
  2:     12050.535: 0.00371847 0.0451786  0.00000  0.00000 0.0575667
0.00000  0.00000 0.00697067 0.0697205 0.0156250 0.0608550 0.0243431
0.0994355
  3:     12037.682: 0.00303460 0.0445012  0.00000  0.00000 0.0568530
0.00000  0.00000 0.00636016 0.0696959 0.0156250 0.0608550 0.0243419
0.0988824
  4:     12012.684: 0.00164710 0.0431313  0.00000  0.00000 0.0554032
0.00000  0.00000 0.00515500 0.0696451 0.0156250 0.0608550 0.0243395
0.0978328
  5:     12003.017: 0.00107848 0.0425739  0.00000  0.00000 0.0548073
0.00000  0.00000 0.00469592 0.0696233 0.0156250 0.0608550 0.0243386
0.0974616
  6:     11984.372:  0.00000 0.0414397  0.00000  0.00000 0.0535899
0.00000  0.00000 0.00378996 0.0695782 0.0156250 0.0608550 0.0243370
0.0967449
  7:     11978.154:  0.00000 0.0409106  0.00000  0.00000 0.0530158
0.00000  0.00000 0.00340746 0.0695560 0.0156250 0.0608550 0.0243363
0.0964537
  8:    -0.0000000:  0.00000      nan  0.00000  0.00000      nan
0.00000  0.00000      nan      nan      nan      nan      nan      nan

Regards,

Jean

2015-05-06 17:43 GMT-07:00 Jean Marchal <jean.d.marchal at gmail.com>:
> 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