From: Douglas Bates <dmbates_at_gmail.com>

Date: Tue 21 Jun 2005 - 23:18:57 EST

> res3 <- nls(y ~ 1/(a*x*x+b*x+c), dd, start = coef(res1), trace = TRUE)

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

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
*