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]]
______________________________________________
[email protected] 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.
