Thanks for the advice! I will continue to monitor the optimizer behaviour. Jean
2015-05-07 17:03 GMT-07:00 William Dunlap <wdun...@tibco.com>: > 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.marc...@gmail.com> > wrote: >> >> Yes, indeed! Problem solved! >> >> Thanks a lot! >> >> Jean >> >> 2015-05-07 14:06 GMT-07:00 William Dunlap <wdun...@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.marc...@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.marc...@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@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. >> > >> > > > ______________________________________________ R-help@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.