R help - Jun 2005 - nls(): Levenberg-Marquardt, Gauss-Newton, plinear - PI curve fitting

Hello,

i have a problem with the function nls().

This are my data in "k":
        V1    V2
 [1,]    0 0.367
 [2,]   85 0.296
 [3,]  122 0.260
 [4,]  192 0.244
 [5,]  275 0.175
 [6,]  421 0.140
 [7,]  603 0.093
 [8,]  831 0.068
 [9,] 1140 0.043
 
With the nls()-function i want to fit following formula whereas a,b, and c
are variables: y~1/(a*x^2+b*x+c)

With the standardalgorithm "Newton-Gauss" the fitted curve contain an
peak
near the second x,y-point.
This peak is not correct for my purpose. The fitted curve should descend
from the maximum y to the minimum y given in my data.

The algorithm "plinear" give me following error:


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

   phi(k[,1],k[,2])

   Error in qr.solve(QR.B, cc) : singular matrix `a' in solve


I have found in the mailinglist
"https://stat.ethz.ch/pipermail/r-help/2001-July/012196.html" that is
if t
he data are artificial. But the data are from my measurment.

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

What are the right way or algorithm for me to solve this problem and what
means this error with alg="plinear"?

Thanks in advance.

-- 
Weitersagen: GMX DSL-Flatrates mit Tempo-Garantie!
Ab 4,99 Euro/Monat: http://www.gmx.net/de/go/dsl

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

On 6/21/05, Christfried Kunath <mailpuls at gmx.net>
wrote:
> Hello,
> 
> i have a problem with the function nls().
> 
> This are my data in "k":
>        V1    V2
>  [1,]    0 0.367
>  [2,]   85 0.296
>  [3,]  122 0.260
>  [4,]  192 0.244
>  [5,]  275 0.175
>  [6,]  421 0.140
>  [7,]  603 0.093
>  [8,]  831 0.068
>  [9,] 1140 0.043
> 
> With the nls()-function i want to fit following formula whereas a,b, and c
> are variables: y~1/(a*x^2+b*x+c)
> 
> With the standardalgorithm "Newton-Gauss" the fitted curve
contain an peak
> near the second x,y-point.
> This peak is not correct for my purpose. The fitted curve should descend
> from the maximum y to the minimum y given in my data.
> 
> The algorithm "plinear" give me following error:
> 
> 
>   phi function(x,y) {
>
k.nls<-nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.0005,b=0.02,c=1.5),alg="plinear")
>       coef(k.nls)
>   }
> 
>   phi(k[,1],k[,2])
> 
>   Error in qr.solve(QR.B, cc) : singular matrix `a' in solve
> 
> 
> I have found in the mailinglist
> "https://stat.ethz.ch/pipermail/r-help/2001-July/012196.html"
that is if t
> he data are artificial. But the data are from my measurment.
> 
> The commercial software "Origin V.6.1" solved this problem with
the
> Levenberg-Marquardt algorithm how i want.
> The reference results are: a = 9.6899E-6, b = 0.00689, c = 2.72982
> 
> What are the right way or algorithm for me to solve this problem and what
> means this error with alg="plinear"?
> 
> Thanks in advance.
This is not a direct answer to your question but log(y) looks nearly linear
in x when plotting them together and log(y) ~ a + b*x or
y ~ a*exp(b*x) will always be monotonic.  Also, this model uses only 2
rather than 3 parameters.