[ADMB Users] NLS example
dave fournier
davef at otter-rsch.com
Fri Feb 22 15:02:11 PST 2013
On 13-02-22 02:27 PM, Ben Bolker wrote:
Its OK Bates has worked it up.
Package source: NISTnls_0.9-13.tar.gz
Here is NLS on the example using the start1 initial values.
> ## starting values for this attempt are ridiculous
> Try(fm1 <- nls(y ~ b1*(x**2+x*b2) / (x**2+x*b3+b4),
+ data = MGH09, trace = TRUE,
+ start = c(b1 = 25, b2 = 39, b3 = 41.5, b4 = 39)))
897.5454 : 25.0 39.0 41.5 39.0
Error in nls(y ~ b1 * (x^2 + x * b2)/(x^2 + x * b3 + b4), data = MGH09, :
singular gradient
and what looks like a different method with the same "ridiculous"
start values
> Try(fm1a <- nls(y ~ b1*(x**2+x*b2) / (x**2+x*b3+b4),
+ data = MGH09, trace = TRUE, alg = "port",
+ start = c(b1 = 25, b2 = 39, b3 = 41.5, b4 = 39)))
0: 448.77269: 25.0000 39.0000 41.5000 39.0000
1: 212.33071: 22.6443 34.8018 45.8794 52.1013
2: 42.041566: 17.2519 25.1976 54.1675 72.0972
3: 2.9195421: 8.94201 14.5357 63.3064 71.4272
4: 0.28575239: 4.69481 19.4082 111.383 113.497
5: 0.14062061: -0.429465 28.6603 125.283 86.6176
6: 0.036909277: 0.784218 16.1017 111.671 203.614
7: 0.0046257983: 0.822336 48.7773 132.310 138.452
8: 0.0020010827: 0.697697 50.8634 141.256 89.4563
9: 0.0015257062: 0.681357 51.1935 142.473 78.3371
10: 0.0010318835: 0.626195 51.8574 143.085 55.2472
11: 0.0010308751: 0.520297 59.4632 141.352 54.8321
12: 0.00097595114: 0.498139 63.2225 140.208 54.6005
13: 0.00096814065: 0.469908 66.6098 139.030 54.3704
14: 0.00096715405: 0.443320 69.7875 137.940 51.9975
15: 0.00095746711: 0.435922 71.3385 137.110 54.1252
16: 0.00095206016: 0.389412 77.6974 134.228 53.0944
17: 0.00093867240: 0.353664 83.8794 130.768 52.0463
18: 0.00092862676: 0.321276 89.8306 126.713 50.6691
19: 0.00091932918: 0.294934 94.9907 122.444 49.1618
20: 0.00091146330: 0.267921 100.324 117.185 47.2623
21: 0.00090298606: 0.244392 105.066 111.521 45.1622
22: 0.00089488572: 0.222512 109.359 105.341 42.8404
23: 0.00089075351: 0.201832 113.046 99.0245 39.0668
24: 0.00087839869: 0.184714 116.278 92.1293 38.1076
25: 0.00086831317: 0.165610 119.225 84.6466 34.8404
26: 0.00086781694: 0.166118 119.225 84.6466 34.8403
27: 0.00086780090: 0.166048 119.226 84.6450 34.8423
28: 0.00086016427: 0.153401 120.962 79.1525 32.8073
29: 0.00085149196: 0.104377 127.007 56.7721 24.1572
30: 0.00073402259: 0.0634784 130.018 33.4309 15.0217
31: 0.00073402259: 0.0634784 130.018 33.4309 15.0217
Error in nls(y ~ b1 * (x^2 + x * b2)/(x^2 + x * b3 + b4), data = MGH09, :
Convergence failure: false convergence (8)
>
Here is ADMB starting from the same "ridiculous" values.
./nlstest -crit 1.e-10 -iprint 50
Initial statistics: 4 variables; iteration 0; function evaluation 0; phase 1
Function value 8.9754538e+02; maximum gradient component mag 7.2704e+01
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 25.00000 7.27040e+01 | 2 39.00000 4.39164e+01 | 3 41.50000
-2.70490e+01
4 39.00000 -1.58952e+01 |
Intermediate statistics: 4 variables; iteration 50; function evaluation
75; phase 1
Function value 9.3002698e-04; maximum gradient component mag 1.3050e-05
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.02942 1.30501e-05 | 2 69.08438 2.10490e-07 | 3 6.20122
5.07298e-08
4 4.22605 7.28610e-08 |
4 variables; iteration 100; function evaluation 136; phase 1
Function value 9.1954376e-04; maximum gradient component mag -9.7628e-03
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.04587 -9.76281e-03 | 2 38.98527 -1.02480e-05 | 3 5.41866
3.83024e-05
4 3.74304 3.83058e-05 |
4 variables; iteration 150; function evaluation 198; phase 1
Function value 8.0592608e-04; maximum gradient component mag -1.6553e-02
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.12474 -1.65531e-02 | 2 5.60738 -2.51737e-04 | 3 1.98041
3.38207e-04
4 1.54576 3.89515e-04 |
- final statistics:
4 variables; iteration 180; function evaluation 236
Function value 3.0751e-04; maximum gradient component mag -9.1546e-11
Exit code = 1; converg criter 1.0000e-10
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.19281 -9.15455e-11 | 2 0.19128 -9.22970e-12 | 3 0.12306
1.08547e-11
4 0.13606 7.38928e-13 |
Estimating row 1 out of 4 for hessian
Estimating row 2 out of 4 for hessian
Estimating row 3 out of 4 for hessian
Estimating row 4 out of 4 for hessian
NLS does converge if you give it really good starting values.
ADMB will converge if you just use the default values of 0
Initial statistics: 4 variables; iteration 0; function evaluation 0; phase 1
Function value 1.4841318e-01; maximum gradient component mag -2.0624e+00
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.00000 -2.06240e+00 | 2 0.00000 0.00000e+00 | 3 0.00000
0.00000e+00
4 0.00000 0.00000e+00 |
ic > imax in fminim is answer attained ?
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.19281 2.63254e-10 | 2 0.19128 8.49052e-11 | 3 0.12306
-5.59366e-11
4 0.13606 -1.28600e-10 |
Function minimizer not making progress ... is minimum attained?
Minimprove criterion = 0.0000e+00
- final statistics:
4 variables; iteration 47; function evaluation 85
Function value 3.0751e-04; maximum gradient component mag 2.6325e-10
Exit code = 1; converg criter 1.0000e-10
Var Value Gradient |Var Value Gradient |Var Value Gradient
1 0.19281 2.63254e-10 | 2 0.19128 8.49052e-11 | 3 0.12306
-5.59366e-11
4 0.13606 -1.28600e-10 |
Estimating row 1 out of 4 for hessian
Estimating row 2 out of 4 for hessian
Estimating row 3 out of 4 for hessian
Estimating row 4 out of 4 for hessian
NLS crashes with the usual ...
> Try(fm2 <- nls(y ~ b1*(x**2+x*b2) / (x**2+x*b3+b4),
+ data = MGH09, trace = TRUE,
+ start = c(b1 = 0, b2 = 0, b3 = 0, b4 = 0)))
Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
What a useless piece of software!
> I *think* (but am not sure) that this is a case where R will do OK --
> it is a case where apparently Excel does something really silly and fits
> by using linear regression on log-transformed data.
>
> I'll work it up if I get a chance.
>
> Ben
>
>
> On 13-02-22 04:41 PM, dave fournier wrote:
>> I'm sure R users would be interested in this!
>>
>>
>>
>> This example runs beautifully with
>>
>> -crit 1.e-10
>>
>> It is part of some excel bashing on the R list in
>>
>> https://stat.ethz.ch/pipermail/r-help/2013-February/347920.html
>>
>> there is na example description
>>
>> http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
>>
>> Wonder how nls does in R
>>
>> tpl,dat and pin files are attached.
>>
>>
>>> On 13-02-22 02:40 PM, dave fournier wrote:
>>>> On 13-02-22 11:33 AM, Ben Bolker wrote:
>>>>
>>>> I agree that backwards compatibility is important.
>>>> We could go back to having a robust flag which would
>>>> activate the options. I think a lot of the existing trouble with this
>>>> data
>>>> set is with
>>>> the bounds on tmpL and u.
>>>> That was probably why sometimes -noinit would crash
>>>> There will always be a new data set
>>>> that will break something.
>>> Of course. I just don't want to break *old* data sets, or at least I
>>> want to allow people to be able to un-break things if the default
>>> settings break them ...
>>>
>>>>> On 13-02-22 02:28 PM, dave fournier wrote:
>>>>>> On 13-02-22 07:26 AM, Ben Bolker wrote:
>>>>>>
>>>>>> Then I reran the model and the variances got stuck near 0.
>>>>>> a small penalty for a while fixes that. So I think a model
>>>>>> something like this would be good.
>>>>>>
>>>>>>
>>>>>>> Hans, can you put this version of the glmmADMB TPL file up on
>>>>>>> the SVN
>>>>>>> repository so the buildbot will build it? It includes a
>>>>>>> beta-binomial
>>>>>>> likelihood option and Dave's latest NB tweak ...
>>>>>>> (I probably do have write access to the SVN repository if I
>>>>>>> could
>>>>>>> figure it out, but I only have to do this a few times a year so I'm
>>>>>>> being lazy ...)
>>>>>>>
>>>>>>> thanks
>>>>>>> Ben
>>>>> Thanks, Dave. I will look this over and try to understand the bits
>>>>> and pieces. If possible I would like to try to make the upper bound on
>>>>> the random effects values into a user-settable parameter (controlled by
>>>>> admb.opts/admbControl) ... Now that we have code that people are
>>>>> using,
>>>>> I worry a lot about making backward-incompatible changes (i.e. this
>>>>> code
>>>>> makes some models run OK, but breaks others that are out there in the
>>>>> wild somewhere).
>>>>>
>>>>> cheers
>>>>> Ben
>>>>>
More information about the Users
mailing list