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

From: Berwin A Turlach <berwin_at_maths.uwa.edu.au>
Date: Tue 21 Jun 2005 - 20:31:04 EST

G'day Chris,

>>>>> "CK" == Christfried Kunath <mailpuls@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@maths.uwa.edu.au Australia http://www.maths.uwa.edu.au/~berwin

R-help@stat.math.ethz.ch mailing list