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

[Gabor Grothendieck]

On 6/21/05, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
> 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.
> 

One other comment.  If you do want to use your model try fitting 1/y
first to get your starting value since that model has a unique solution:

> res1 <- nls(1/y ~ a*x^2 + b*x + c, start = list(a=0,b=0,c=0))
> res2 <- nls(y ~ 1/(a*x^2 + b*x + c), start = as.list(coef(res1)))
> res2
Nonlinear regression model
  model:  y ~ 1/(a * x^2 + b * x + c) 
   data:  parent.frame() 
           a            b            c 
9.690187e-06 6.885577e-03 2.729825e+00 
 residual sum-of-squares:  0.000547369