Dear Nilsson,
Thank you very much for your reply which gave me much instruction!
I will have a try on the suggested way from you.
In fact, to get the value of MA is the most important thing for me
in the fitting by this model. And I have thought about the linear model
fitting at first, but you know, when the the functional form be changed
as:
y = (a + b*c^2 + d*c^4) + (-2*b*c - 4*d*c^3)*x +
(d + b*6*c^2)*x^2 - (4*d*c)*x^3 + d*x^4 + epsilon
four coefficient (a,b,c,d) turns to five in this linear model (for
constant, x, x^2, x^3, x^4), that means there should be some
relationship among the five coefficients, which cannot be reflected in
the linear model. So I cannot get the accurate value for the parameters.
For your reference:
As what you pointed out, the initial value for a,b,c,MA is critical
for nls. In recent days, I tried recruiting the model, y~a*x^4+b*x^2+c,
to get the initial value of a0, b0, and c0, and it works, but cannot
avoid the convergence problem totally (the probability went down).
Thank you very much again!
I really appreciate your help!
Best regards,
Warren
Nilsson Fredrik X wrote:
> Dear Warren,
>
> I had similar problems, this is (roughly) how I solved it translated to your
> problem:
> x <- c(-0.06,-0.04,-0.025,-0.015,-0.005,0.005,0.015,0.025,0.04,0.06)
>
> y <- c(1866760,1457870,1314960,1250560,1184850,1144920,
> 1158850,1199910,1263850,1452520)
>
> dafa<-data.frame(x,y)
> plot(x,y)
>
> foqufu<-function(parm)
> {
> const<-parm[1]
> A<-parm[2]
> B<-parm[3]
> MA<-parm[4]
>
> d<-y-(const + A*(x-MA)^4 + B*(x-MA)^2)
> d<-d%*%d
> return(d)
> }
>
> # This is your initial guess
> aaa<-c(10000000, 100000000, -1000000, 0)
> # As already pointed out, you have a bad initial guess.
>
>
> apa3<-optim(aaa, foqufu, control=list(maxit=20000, reltol=1e-16))
> apa3<-optim(apa3$par, foqufu, control=list(maxit=20000, reltol=1e-16))
> apa3<-optim(apa3$par, foqufu, control=list(maxit=20000, reltol=1e-16))
> apa3<-optim(apa3$par, foqufu, control=list(maxit=20000, reltol=1e-16))
> # I do this a couple of times until convergence
> apa3<-optim(apa3$par, foqufu, control=list(maxit=20000, reltol=1e-10))
>
> fitOup<- nls(y ~ constant + A*(x-MA)^4 + B*(x-MA)^2, data=dafa,
> start=list(constant= 4.911826e+06, A=2.625077e+08, B=-6.278897e+07,
> MA=3.538064e-01),
> control=nls.control(maxiter=100, minFactor=1/4096), trace=TRUE)
>
>
> Note that due to your bad initial guess on parameters you end up in a local
> not global minimum which is the trouble with nonlinear optimization. A way
> out would perhaps be to randomize your initial guesses.
>
> Andy's way of getting the initial values is better in that respect. NOTE that
> he searches for the quadratic first (downplaying the importance of the
> quartic term) But technically I don't understand the need to scale your y's,
> in this case you get the same.
>
> Another way of finding some intial values is to solve (const + B*MA^2 =
> intercept, -2*B*MA = coefficient of x, B = coefficient of x^2 in the J2.lm
> fit below; assuming that J2.lm (below has been fit):
> B<-as.numeric(coef(J2.lm)[3])
> A<-0
> MA<-as.numeric(coef(J2.lm)[2])/2/B
> const<- as.numeric(coef(J2.lm)[1])-B*MA^2
>
> bbb<-c(const, 0, B, MA)
> apa3<-optim(bbb, foqufu, control=list(maxit=20000, reltol=1e-16))
> #iterating a couple of times
> apa3<-optim(apa3$par, foqufu, control=list(maxit=20000, reltol=1e-16))
> #gives the same solution as Andy's
> ).
> #variation of Andy's solution (no scaling of y)
> J2.lm<-lm(y~1 + x + I(x^2), data=dafa)
> summary(J2.lm)
>
> plot(xx, predict(J2.lm, ydf))
> xx[which.min(predict(J2.lm, ydf))]
> Ja2.lm<-lm(y~1 + I((x-0.01)^2), data=dafa)
> summary(Ja2.lm)
> plot(x,y)
> points(x, predict(Ja2.lm, data=data.frame(x=x-0.01, y=y)), col="red")
>
> apa <-optim(c(1146530, 0,139144223, 0.01), foqufu, control=list(maxit=20000,
> reltol=1e-10))
>
> fitOup3b<- nls(y ~ constant + A*(x-MA)^4 + B*(x-MA)^2, data=dafa,
> start=list(constant= apa$par[1], A=apa$par[2], B= apa$par[3],
> MA=apa$par[4]),
> control=nls.control(maxiter=100, minFactor=1/4096), trace=TRUE)
> # this is exactly the same solution as Andy's, except that the residuals are
> # 1e12 larger and A,B, const are 1e6 larger (due to the scaling).
>
>
> BUT, and this is a big BUT, why on earth do you choose this functional form?
> First, the plot(x,y) suggest to my eyes that the function is not symmetric
> around its centre. Second, if you're not particularly interested in finding a
> value for MA (i.e. it can be considered as a constant) why not use a simple
> linear regression with a polynomial? What I mean here is that if:
>
> y = a + b*(x-c)^2 + d*(x-c)^4 + epsilon
> (and no variability in c) then
> y = (a + b*c^2 + d*c^4) + (-2*b*c - 4*d*c^3)*x +
> (d + b*6*c^2)*x^2 - (4*d*c)*x^3 + d*x^4 + epsilon
>
> So in fact you could do:
> A.lm<-lm(y ~ poly(x,4), data=dafa)
> Summary(A.lm)
> # suggests that fourth order not necessary
> B.lm<-lm(y ~ poly(x,3), data=dafa).
> anova(B.lm, A.lm)
>
> # note that this suggest that your problem has an asymmetry that
> # should be accounted for (unless you have very strong reasons for choosing
> # your particular model formulation).
>
> Third, and this is a more critical question (which I hope that the experts on
> nls/R gurus would comment on how to solve) I think it is problematic to have
> variability in LOCATION such as your shift variable (MA). It seems to
> introduce all sorts of nastiness, which ought to be particularly severe in
> non-linear regression.
>
> Best regards,
> Fredrik
>
> -----Ursprungligt meddelande-----
> Från: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> Skickat: den 5 september 2007 18:24
> Till: Yu (Warren) Wang
> Kopia: [EMAIL PROTECTED]; [EMAIL PROTECTED]
> Ämne: Re: [R] question about non-linear least squares in R
>
> Here is one way of getting a reasonable fit:
>
> 1. Scale your y's by dividing all values by 1e6.
> 2. Plot x vs. y. The plot looks like a quadratic function.
> 3. Fit a quadratic const. + B*x^2 - this a linear regression problem so
> use lm.
> 4. Plot the predictions.
> 5. Eyball the necessary shift - MA is around 0.01. Refit const. +
> B*(x-.01)^2. Should get const.=1.147 and B=139.144
> 6. Use start=list(const.= 1.147, A=0, B=1.147, MA=.01). nls should
> converge in 4 iterations.
>
> In general, good starting points may be crucial to nls convergence.
> Scaling the y's to reasonable values also helps.
>
> Hope this helps,
>
> Andy
>
> __________________________________
> Andy Jaworski
> 518-1-01
> Process Laboratory
> 3M Corporate Research Laboratory
> -----
> E-mail: [EMAIL PROTECTED]
> Tel: (651) 733-6092
> Fax: (651) 736-3122
>
>
>
> "Yu (Warren)
> Wang"
> <[EMAIL PROTECTED]> To
> Sent by: "[EMAIL PROTECTED]"
> [EMAIL PROTECTED] <[EMAIL PROTECTED]>
> at.math.ethz.ch cc
>
> Subject
> 09/05/2007 02:51 [R] question about non-linear least
> AM squares in R
>
>
>
>
>
>
>
>
>
>
> Hi, everyone,
> My question is: It's not every time that you can get a converged
> result from the nls function. Is there any solution for me to get a
> reasonable result? For example:
>
> x <- c(-0.06,-0.04,-0.025,-0.015,-0.005,0.005,0.015,0.025,0.04,0.06)
>
> y <-
> c(1866760,1457870,1314960,1250560,1184850,1144920,1158850,1199910,1263850,1452520)
>
>
> fitOup<- nls(y ~ constant + A*(x-MA)^4 + B*(x-MA)^2,
> start=list(constant=10000000, A=100000000, B=-1000000, MA=0),
> control=nls.control(maxiter=100, minFactor=1/4096), trace=TRUE)
>
>
>
> For this one, I cannot get the converged result, how can I reach it? To
> use another funtion or to modify some settings for nls?
>
> Thank you very much!
>
> Yours,
>
> Warren
>
> ______________________________________________
> [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.
>
>
>
>
______________________________________________
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.
