Re: [R] nls.control

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Wed 06 Apr 2005 - 22:21:49 EST

>>>>> "Anthony" == Anthony Landrevie <antho_l@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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Apr 06 22:25:28 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:31:02 EST