[R] Upper limit in nlsLM not working as expected

[Martin Hehn]

Thanks for the reply Berend,

I am aware that I do not have to define the limits to Inf or -Inf, I just did this
to make sure all other variables (besides 'w') in 'upper' have no limits.

I can tolerate no upper limit as my peak is fitted well, but my actual data is a time series where the peak disappears from time to time. Once disappeared, the algorithm is still fitting to the noise of the experimental dataset. 

I thought that it would be handy to suppress this by simply specifying a limit for how wide I allow the peak to be as it becomes very broad when fitting to noise. Any ideas on how to avoid fitting to the noise?

Regards 
________________________________________
From: Berend Hasselman [bhh at xs4all.nl]
Sent: 18 October 2012 17:45
To: Martin Hehn
Cc: r-help at R-project.org
Subject: Re: [R] Upper limit in nlsLM not working as expected

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.