R help - May 2015 - nlminb supplying NaN parameters to objective function

Your nLL function returns 1e+308 in near-boundary cases.  Since 1e+308 is so
close to machine infinity, it is easy to get into Inf-Inf (=NaN) or Inf/Inf
(=NaN)
situations when working with it.  Try making that limiting value something
smaller,
like 1e+30, and you may have better luck.

Bill Dunlap
TIBCO Software
wdunlap tibco.com

On Thu, May 7, 2015 at 1:14 PM, Jean Marchal <jean.d.marchal at gmail.com>
wrote:
> 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
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]

Yes, indeed! Problem solved!

Thanks a lot!

Jean

2015-05-07 14:06 GMT-07:00 William Dunlap <wdunlap at
tibco.com>:
> Your nLL function returns 1e+308 in near-boundary cases.  Since 1e+308 is
so
> close to machine infinity, it is easy to get into Inf-Inf (=NaN) or Inf/Inf
> (=NaN)
> situations when working with it.  Try making that limiting value something
> smaller,
> like 1e+30, and you may have better luck.
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
> On Thu, May 7, 2015 at 1:14 PM, Jean Marchal <jean.d.marchal at
gmail.com>
> wrote:
>>
>> 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
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>

Your immediate problem may be solved, but the exact value of that limiting
value
affects the parameter estimates a fair bit.  I have not really looked at
your function,
but the ledge around it puts a kink (discontinuous first derivative) into
it, which can
mess up optimizers.

Bill Dunlap
TIBCO Software
wdunlap tibco.com

On Thu, May 7, 2015 at 4:46 PM, Jean Marchal <jean.d.marchal at gmail.com>
wrote:
> Yes, indeed! Problem solved!
>
> Thanks a lot!
>
> Jean
>
> 2015-05-07 14:06 GMT-07:00 William Dunlap <wdunlap at tibco.com>:
> > Your nLL function returns 1e+308 in near-boundary cases.  Since 1e+308
> is so
> > close to machine infinity, it is easy to get into Inf-Inf (=NaN) or
> Inf/Inf
> > (=NaN)
> > situations when working with it.  Try making that limiting value
> something
> > smaller,
> > like 1e+30, and you may have better luck.
> >
> > Bill Dunlap
> > TIBCO Software
> > wdunlap tibco.com
> >
> > On Thu, May 7, 2015 at 1:14 PM, Jean Marchal <jean.d.marchal at
gmail.com>
> > wrote:
> >>
> >> 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
> >>
> >> ______________________________________________
> >> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more,
see
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >
> >
>
	[[alternative HTML version deleted]]