[R] Upper limit in nlsLM not working as expected

[nashjc at uottawa.ca]

I think this might be what you want.  Kate Mullen and I have been in
correspondence over some edge cases where minpack.LM may not handle bounds
appropriately. However, though nlmrt seems to do the job here, readers
should note that R benefits hugely if we maintain some friendly
competition (and collaboration) to both offer methods with different
performance profiles (e.g., handling different classes/sizes of problems
better) and in
helping to root out bugs. So I want minpack.LM to be there for comparison
and use when itdoes a better job. It should, for example, be much faster.
nlmrt is designed to be aggressive and also modifiable, so that the
algorithm can be explored.

John Nash


library(nlmrt)
#Fit function to data
fitn<-nlxb(y ~ y0 + 2*(A/pi)*w/(4*(x-xc)^2+w^2),
            control=nls.lm.control(
              factor=100,
              maxiter=1024,
              ftol = .Machine$double.eps,
              ptol = .Machine$double.eps
            ),
            data=data,
            na.action=na.exclude,
            start=st2,
            algorith='LM',
            lower=c(-0.0001,-1e-8,0.05,0.2),
            #upper=c(1e-6,0.003,0.08,0.35),
            upper=c(Inf,Inf,0.07,Inf),
            trace=F
)
# > str(fitn)
# List of 6
#  $ resid   : num [1:151] 0.00265 0.00124 -0.00222 -0.00239 0.00148 ...
#  $ jacobian: num [1:151, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
#   ..- attr(*, "dimnames")=List of 2
#   .. ..$ : NULL
#   .. ..$ : chr [1:4] "y0" "A" "w" "xc"
#  $ feval   : num 24
#  $ jeval   : num 16
#  $ coeffs  : num [1:4] 0.00014 0.0099 0.05939 0.22975
#  $ ssquares: num 0.000857




On 10/19/2012 06:00 AM, r-help-request at r-project.org wrote:
> Message: 11
> Date: Thu, 18 Oct 2012 12:16:17 +0000
> From: Martin Hehn <Martin.Hehn at postgrad.manchester.ac.uk>
> To: "r-help at R-project.org" <r-help at R-project.org>
> Subject: [R] Upper limit in nlsLM not working as expected
> Message-ID:
> 	<17F41A7342EEB6409BDA31C3ACDBAC1005AD55 at MBXP07.ds.man.ac.uk>
> Content-Type: text/plain
>
> Dear all,
>
> I am using the nlsLM function to fit a Lorentzian function to my
experimental data.
> The LM algorithm should allow to specify limits, but the upper limit
appears not to work as expected in my code.
> The parameter 'w', which is peak width at half maximuim always hits the
upper limit if the limit is specified. I would expect the value to be
in-between the upper and lower limit with maybe hitting either limit
occasionally.
>
> The code below calculates a lorentzian curve with some noise added that
looks like a typical experimental spectrum. Then a fit is made to this
data. Note that if 'w' in the upper limit is set to e.g 0.6, 0.7, or 0.8
the fit always hits this limit. If 'w'=Inf, the fit calculates to
something around 0.06, which is correct.
>
> Why does the fit go wrong if 'w' is not Inf?
>
> library(minpack.lm)
> #Create x axis
> x<-seq(from=0.1,to=0.6,by=0.5/150)
> #Simulate a function with noise
> fu<-function(y0,A,w,xc,x){
>   eq<-y0 + 2*A/pi*w/(4*(x-xc)^2+w^2)
>   eq<-jitter(eq,factor=200)
>   }
> #Evaluate function aka Measured data
> y<-fu(0,0.01,0.06,0.23,x)
> data<-as.data.frame(cbind(x,y))
>
> #Start values for fitting
> st2<-data.frame(
>   y0=0,
>   A=0.0001,
>   w=0.055,
>   xc=0.28
> )
> #Fit function to data
> fit<-nlsLM(y ~ y0 + 2*(A/pi)*w/(4*(x-xc)^2+w^2),
>             control=nls.lm.control(
>               factor=100,
>               maxiter=1024,
>               ftol = .Machine$double.eps,
>               ptol = .Machine$double.eps
>             ),
>             data=data,
>             na.action=na.exclude,
>             start=st2,
>             algorith='LM',
>             lower=c(-0.0001,-1e-8,0.05,0.2),
>             #upper=c(1e-6,0.003,0.08,0.35),
>             upper=c(Inf,Inf,0.07,Inf),
>             trace=F
> )
> #Predict fitting values
> fity<-predict(fit,data$x)
>
> plot(data$x,data$y)
> lines(data$x,fity,col=2)
> text(0.4,0.08,coef(fit)['y0'])
> text(0.4,0.07,coef(fit)['A'])
> text(0.4,0.06,coef(fit)['w'])
> text(0.4,0.05,coef(fit)['xc'])
>
> Best regard
> Martin