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