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

From: Manuel Morales <Manuel.A.Morales_at_williams.edu>
Date: Wed 22 Jun 2005 - 00:22:56 EST

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

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



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 Wed Jun 22 00:32:16 2005

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