[R] nls.control

Martin Maechler maechler at stat.math.ethz.ch
Wed Apr 6 14:21:49 CEST 2005


>>>>> "Anthony" == Anthony Landrevie <antho_l at yahoo.com>
>>>>>     on Wed, 6 Apr 2005 02:54:50 -0700 (PDT) writes:

    Anthony> Hello everyone, I'm trying to test the accurracy of
    Anthony> R on the Eckerle4 dataset from NIST 

Do you know that there's an R package 'NISTnls' (on CRAN)
exactly for this purpose?
After installing the package,

  library(NISTnls)
  example(Eckerle4)

gives

  Eckrl4> data(Eckerle4)

  Eckrl4> plot(y ~ x, data = Eckerle4)

  Eckrl4> fm2 <- nls(y ~ (b1/b2) * exp(-0.5 * ((x - b3)/b2)^2), 
      data = Eckerle4, trace = TRUE, start = c(b1 = 1.5, b2 = 5, 
	  b3 = 450))
  0.05668291 :    1.5   5.0 450.0 
  0.00722609 :    1.563149   4.374689 451.974368 
  0.001525831 :    1.551040   4.091636 451.488425 
  0.001463731 :    1.554819   4.091467 451.541251 
  0.001463589 :    1.554395   4.088899 451.541108 
  0.001463589 :    1.554384   4.088839 451.541216 
  0.001463589 :    1.554383   4.088832 451.541218 

  Eckrl4> fm4 <- nls(y ~ (1/b2) * exp(-0.5 * ((x - b3)/b2)^2), 
      data = Eckerle4, trace = TRUE, start = c(b2 = 5, b3 = 450), 
      algorithm = "plinear")
  0.05086068 :   5.00000 450.00000   1.65696 
  0.004539377 :   4.471095 451.669974   1.621837 
  0.001478679 :   4.085508 451.514686   1.553734 
  0.001463615 :   4.089948 451.541333   1.554595 
  0.001463589 :   4.088856 451.541172   1.554387 
  0.001463589 :   4.088835 451.541217   1.554383 
  0.001463589 :   4.088832 451.541218   1.554383 
  > 
 --------------

where the "fm2 <- " case looks pretty much like your example below

    Anthony> and I don't understand how the control option of
    Anthony> the nls function works.  I tought nls(...) was
    Anthony> equivalent to nls(...control=nls.control()) i.e
    Anthony> nls.control() was the default value of control, but
    Anthony> here is the error I get :
 
    >> n2=nls(V1~(b1/b2) *
    >> exp(-0.5*((V2-b3)/b2)^2),data=ecker,start=list(b1=1.5,b2=5,b3=450,control=nls.control()))
    Anthony> Error in nlsModel(formula, mf, start) : singular
    Anthony> gradient matrix at initial parameter estimates

we cannot know, since we don't have your "ecker".

For me, with 'Eckerle4' from the "NISTnls" package,

  fm2. <- nls(y ~ (b1/b2) * exp(-0.5 * ((x - b3)/b2)^2), data = Eckerle4, trace = TRUE, start = c(b1 = 1.5, b2 = 5, b3 = 450), control=nls.control())

  0.05668291 :    1.5   5.0 450.0 
  0.00722609 :    1.563149   4.374689 451.974368 
  0.001525831 :    1.551040   4.091636 451.488425 
  0.001463731 :    1.554819   4.091467 451.541251 
  0.001463589 :    1.554395   4.088899 451.541108 
  0.001463589 :    1.554384   4.088839 451.541216 
  0.001463589 :    1.554383   4.088832 451.541218 

I get exactly the same when I have added  
  " , control=nls.control() "
to the original call.  So I wonder if you didn't accidentally
change something more than just adding that.

    Anthony> while I get no error without setting the control
    Anthony> option with the same other parameters.
 
    Anthony> I see that R didn't manage 

that's a pretty tough statement (and really wrong). I assume it should
mean 
     "nls() didn't solve ...., at least not with default
      arguments specified"

and I think you are right: completely wrong starting values
don't always "work" for nls()

    Anthony> to solve the Eckerle4 regression problem from start one

"start one" is the infamous non-sense of obviously completely
wrongly specified starting values, right?

    Anthony> while Splus can do it with the nlregb option.  Is there something
    Anthony> equivalent for R now?
 
not "equivalent" probably, but yes, there are several
alternatives for minimization/optimization, in 
"base+recommended R", and more in other packages.


    Anthony> Otherwise, I found that R 2.0.1 was performing
    Anthony> better than SAS 9.1 on the NIST Datasets in
    Anthony> general.
 
where  "R 2.0.1"  means using nls() in R 2.0.1, right?
Thanks for letting us know.

    Anthony> Best regards,
 
    Anthony> Anthony Landrevie




More information about the R-help mailing list