[R] nls(): Levenberg-Marquardt, Gauss-Newton, plinear - PI curve fitting

[Berwin A Turlach]

G'day Chris,

>>>>> "CK" == Christfried Kunath <mailpuls at gmx.net> writes:

    CK> With the nls()-function i want to fit following formula
    CK> whereas a,b, and c are variables: y~1/(a*x^2+b*x+c)

    CK> [...]

    CK> The algorithm "plinear" give me following error:
The algorithm "plinear" is inappropriate for your data and your model
since none of the parameters are linear.  You are actually trying to
fit the model
                y ~ d/(a*x^2+b*x+c)

where `d' would be the linear parameter.

    CK> phi function(x,y) {
    CK> k.nls<-nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.0005,b=0.02,c=1.5),alg="plinear")
    CK> coef(k.nls) }

    CK> [...]

    CK> The commercial software "Origin V.6.1" solved this problem
    CK> with the Levenberg-Marquardt algorithm how i want.  The
    CK> reference results are: a = 9.6899E-6, b = 0.00689, c = 2.72982

    CK> What are the right way or algorithm for me to solve this
    CK> problem and what means this error with alg="plinear"?
The error means that at some point along the way a matrix was
calculated that needed to be inverted but was for all practical
purposes singular.  This can happen in numerical optimisation
problems, in particular if derivatives have to be calculated
numerically.

How to solve this problem:

1) Don't use the algorithm "plinear" since it is inappropriate for
   your model.

2) You may want to specify the gradient of the function that you are
   minimising to make life easier for nls(), see Venables & Ripley
   (2002, page 215) for an example.

3) You can call nls() directly without specifying the plinear option:
   (I renamed the variables in the data frame to x and y for
   simplicity) 

      > nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.0005,b=0.02,c=1.5),data=k)
      Nonlinear regression model
        model:  y ~ 1/(a * (x^2) + b * x + c) 
         data:  k 
                  a             b             c 
      -7.326117e-05  4.770514e-02  2.490643e+00 
       residual sum-of-squares:  0.1120086

   But the results seem to be highly depended on your starting values:

      > nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.00005,b=0.002,c=2.5),data=k)
      Nonlinear regression model
        model:  y ~ 1/(a * (x^2) + b * x + c) 
         data:  k 
                 a            b            c 
      9.690204e-06 6.885570e-03 2.729825e+00 
       residual sum-of-squares:  0.000547369 

   Which is of some concern.

4) If you really want to fit the above model, you may also consider to
   just use the glm() command and fit it within a generalised linear
   model framework:

      > glm(y ~ I(x^2) + x, data=k, family=gaussian(link="inverse"))
      
      Call:  glm(formula = y ~ I(x^2) + x, family = gaussian(link = "inverse"),      data = k) 
      
      Coefficients:
      (Intercept)       I(x^2)            x  
        2.730e+00    9.690e-06    6.886e-03  
      
      Degrees of Freedom: 8 Total (i.e. Null);  6 Residual
      Null Deviance:	    0.09894 
      Residual Deviance: 0.0005474 	AIC: -53.83 

HTH.

Cheers,

        Berwin

========================== Full address ============================
Berwin A Turlach                      Tel.: +61 (8) 6488 3338 (secr)   
School of Mathematics and Statistics        +61 (8) 6488 3383 (self)      
The University of Western Australia   FAX : +61 (8) 6488 1028
35 Stirling Highway                   
Crawley WA 6009                e-mail: berwin at maths.uwa.edu.au
Australia                        http://www.maths.uwa.edu.au/~berwin