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

[Douglas Bates]

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

In cases where you are having difficulty getting nls to converge
(which is not the case the way that Gabor is fitting the model) it is
best to use trace = TRUE to see what is happening to the parameter
estimates.  You can use the alg = 'plinear' option for this model if
you divide the numerator and denominator of the rational function by
one of the parameters.  For example, if you divide by a then you have
a monic quadratic in the denominator and a parameter equivalent to 1/a
in the numerator.  I chose to divide by c but it could work with any
of the parameters

> res1 <- nls(1/y ~ a*x*x + b*x + c, dd, start=c(a=0, b=0, c=0), trace = TRUE)
1006.817 :  0 0 0 
0.597423 :  9.751120e-06 6.755636e-03 2.760013e+00 
> coef(res1)
           a            b            c 
9.751120e-06 6.755636e-03 2.760013e+00 
> coef(res1)[-3]/coef(res1)[3]
           a            b 
3.532998e-06 2.447683e-03 
> res2 <- nls(y ~ 1/(a*x*x+b*x+1), dd, start = coef(res1)[-3]/coef(res1)[3], alg = 'plinear', trace = TRUE)
0.0005553948 : 3.532998e-06 2.447683e-03 3.643117e-01 
0.0005473694 : 3.549814e-06 2.521742e-03 3.663088e-01 
0.000547369 : 3.549764e-06 2.522343e-03 3.663238e-01 
> res3 <- nls(y ~ 1/(a*x*x+b*x+c), dd, start = coef(res1), trace = TRUE)
0.0005678106 :  9.751120e-06 6.755636e-03 2.760013e+00 
0.0005473708 :  9.690171e-06 6.886612e-03 2.729552e+00 
0.000547369 :  9.690188e-06 6.885577e-03 2.729825e+00 

The solutions are equivalent.

> coef(res3)[-3]/coef(res3)[3]
           a            b 
3.549747e-06 2.522351e-03