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
You have nearly as many parameters as data points which may cause fundamental problem singularity problems but one thing to try, just in case, is to transform it to unconstrained. For example, let A1 = 0.1 + A1x^2 and then substitute A1 with the right hand side so that it becomes a function of the parameter A1x instead of A1 and similarly for the other parameters so that it becomes unconstrained. On 5/21/06, Lorenzo Isella <lorenzo.isella at gmail.com> 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 at 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 >
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 at 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.htmlapart 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).