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

From: Douglas Bates <dmbates_at_gmail.com>
Date: Tue 21 Jun 2005 - 23:18:57 EST

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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Jun 21 23:21:48 2005

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