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

[Douglas Bates]

On 6/21/05, Manuel Morales <Manuel.A.Morales at williams.edu> wrote:
> 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.

One reason that will be slow is because you are doing an unnecessary
loop in your fit.model function.  You can write the sum of squares
function as

function(pars) sum((y - 1/(pars[1]*x*x + pars[2]*x + pars[3]))^2)

Remember that arithmetic operations in the S language are vectorized.