[R] Upper limit in nlsLM not working as expected

[Berend Hasselman]

On 18-10-2012, at 14:16, Martin Hehn wrote:

> 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?
> 

The upper limit for 'w' does not need to be set to Inf. It needs to be set high enough for LM to 'work'.
If you set the upper limit for 'w' to 10 the algorithm  will also estimate 'w' at 0.06..

When you set the upper limit to low values as you have done, you are not giving LM enough room to manoeuvre.
If you set trace to TRUE (please don't use F and T) you'll see that 'w' hits the upper limit and that it stays there.

Berend


> 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
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.