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