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
[[alternative HTML version deleted]]
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.
nashjc at uottawa.ca
2012-Oct-19 12:50 UTC
[R] Upper limit in nlsLM not working as expected
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
Dear John, this appears to be one of the cases where your algorithm is favourable over the nlsLM one. The fitting is now a lot better with ssquared being lower compared to any other fit method I was comparing to. As I am relatively new to R, it still puzzles me that there are so many packages doing what I want but they seem so hard to track down. Best regards, Martin