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

[Manuel Morales]

On Tue, 2005-06-21 at 06:57 -0400, Gabor Grothendieck 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.

If you want to use the original model, you could also use optim() and
specify your own minimization function. The default algorithm is slow
but has worked well for me where the faster algorithms in nls() have
choked. I adapted one of my functions for your data below:

minimize.fn<-function(parms) {

  fit.model<-function(parms) {
    y.pred={}; for (i in 1:9) {
      a=parms[1];b=parms[2];c=parms[3];
      y=1/(a*V1[i]^2+b*V1[i]+c)
      y.pred=append(y.pred,y)}
    list(fitted.values=y.pred)}

  minimize.model<-function(parms) {
    Resids=V2-fit.model(parms)$fitted.values
    ss=sum(Resids^2)
    list(c(ss))}

  fit=optim(c(parms[1],parms[2],parms[3]),minimize.model,
    control=list(trace=1,maxit=10000))
  results=fit.model(c(fit$par[1],fit$par[2],fit$par[3]))
 
  list(parms=fit$par,SS=fit$value,residuals=results$fitted.values-V2,
    fitted.values=results$fitted.values)}

result<-minimize.fn(c(0.0005,.02,1.5))
result<-minimize.fn(result$parms)

This gives the results:

> > result$parms
> [1] 1.184172e-05 4.878992e-03 3.045663e+00
> > result$SS
> [1] 0.005436355

Which is very similar to what you want.

HTH,

Manuel