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

From: Berwin A Turlach <>
Date: Tue 21 Jun 2005 - 20:31:04 EST

G'day Chris,

>>>>> "CK" == Christfried Kunath <> 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) 
      (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 mailing list PLEASE do read the posting guide! Received on Tue Jun 21 20:39:21 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:55 EST