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-requ...@r-project.org wrote:
> Message: 11
> Date: Thu, 18 Oct 2012 12:16:17 +0000
> From: Martin Hehn <martin.h...@postgrad.manchester.ac.uk>
> To: "r-help@R-project.org" <r-help@R-project.org>
> Subject: [R] Upper limit in nlsLM not working as expected
> Message-ID:
> <17f41a7342eeb6409bda31c3acdbac1005a...@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
______________________________________________
R-help@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.
