From: Joerg van den Hoff <j.van_den_hoff_at_fz-rossendorf.de>

Date: Mon 22 May 2006 - 17:57:26 EST

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon May 22 18:02:13 2006

Date: Mon 22 May 2006 - 17:57:26 EST

Lorenzo Isella wrote:

> Dear All,

*> I may look ridiculous, but I am puzzled at the behavior of the nls with
**> a fitting I am currently dealing with.
**> My data are:
**>
**> x N
**> 1 346.4102 145.428256
**> 2 447.2136 169.530634
**> 3 570.0877 144.081627
**> 4 721.1103 106.363316
**> 5 894.4272 130.390552
**> 6 1264.9111 36.727069
**> 7 1788.8544 52.848587
**> 8 2449.4897 25.128742
**> 9 3464.1016 7.531766
**> 10 4472.1360 8.827367
**> 11 6123.7244 6.600603
**> 12 8660.2540 4.083339
**>
**> I would like to fit N as a function of x according to a function
**> depending on 9 parameters (A1,A2,A3,mu1,mu2,mu3,myvar1,myvar2,myvar3),
**> namely
**> N ~
**> (log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1))
**>
**> +log(10)*A2/sqrt(2*pi)/log(myvar2)*exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2))
**>
**> +log(10)*A3/sqrt(2*pi)/log(myvar3)*exp(-((log(x/mu3))^2)/2/log(myvar3)/log(myvar3)))
**>
**> (i.e. N is to be seen as a sum of three "bells" whose parameters I need
**> to determine).
**>
**>
**> So I tried:
**> out<-nls(N ~
**> (log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1))
**>
**> +log(10)*A2/sqrt(2*pi)/log(myvar2)*exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2))
**>
**> +log(10)*A3/sqrt(2*pi)/log(myvar3)*exp(-((log(x/mu3))^2)/2/log(myvar3)/log(myvar3)))
**> ,start=list(A1 = 85,
**> A2=23,A3=4,mu1=430,mu2=1670,mu3=4900,myvar1=1.59,myvar2=1.5,myvar3=1.5 )
**> ,algorithm = "port"
**> ,control=list(maxiter=20000,tol=10000)
**> ,lower=c(A1=0.1,A2=0.1,A3=0.1,mu1=0.1,mu2=0.1,mu3=0.1,myvar1=0.1,myvar2=0.1,myvar3=0.1)
**> )
**>
**> getting the error message:
**> Error in nls(N ~ (log(10) * A1/sqrt(2 * pi)/log(myvar1) *
**> exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1)) + :
**> Convergence failure: singular convergence (7)
**>
**>
**> I tried to adjust tol & maxiter, but unsuccessfully.
**> If I try fitting N with only two "bells", then nls works:
**>
**> out<-nls(N ~
**> (log(10)*A1/sqrt(2*pi)/log(myvar1)*exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1))
**>
**> +log(10)*A2/sqrt(2*pi)/log(myvar2)*exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2))
**> )
**> ,start=list(A1 = 85, A2=23,mu1=430,mu2=1670,myvar1=1.59,myvar2=1.5 )
**> ,algorithm = "port"
**> ,control=list(maxiter=20000,tol=10000)
**> ,lower=c(A1=0.1,A2=0.1,mu1=0.1,mu2=0.1,myvar1=0.1,myvar2=0.1)
**> )
**>
**> out
**> Nonlinear regression model
**> model: N ~ (log(10) * A1/sqrt(2 * pi)/log(myvar1) *
**> exp(-((log(x/mu1))^2)/2/log(myvar1)/log(myvar1)) + log(10) *
**> A2/sqrt(2 * pi)/log(myvar2) *
**> exp(-((log(x/mu2))^2)/2/log(myvar2)/log(myvar2)))
**> data: parent.frame()
**> A1 A2 mu1 mu2 myvar1 myvar2
**> 84.920085 40.889968 409.656404 933.081936 1.811560 2.389215
**> residual sum-of-squares: 2394.876
**>
**> Any idea about how to get nls working with the whole model?
**> I had better luck with the nls.lm package, but it does not allow to
**> introduce any constrain on my fitting parameters.
**> I was also suggested to try other packages like optim to do the same
**> fitting, but I am a bit unsure about how to set up the problem.
**> Any suggestions? BTW, I am working with R Version 2.2.1
**>
**> Lorenzo
**>
**> ______________________________________________
**> R-help@stat.math.ethz.ch mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-help
**> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
*

apart from the fact that fitting 9 parameters to 12 points quite genereally will not yield satisfactory results (at least estimates will have huge uncertainties), total failure in your case seems obviouus if you plot your data: it's not even obvious where the three peaks (means of the gaussians) should be: all below x=2000 or is there one peak at about x=4500 and one of the 'peaks' below x=2000 is spurious? if you can't decide, nls has problems. moreover: how should reliable estimates of the standard deviations (width) of the gaussian result if the peaks essentially consist of exactly one point?

in short: I believe, you either need much more data points or you must put in substantial a priori information (e.g. either means or standard deviations of the gaussians).

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon May 22 18:02:13 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Mon 22 May 2006 - 20:10:17 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*