[R] question about non-linear least squares in R
Yu (Warren) Wang
yu.wang at pdf.com
Wed Sep 12 03:14:23 CEST 2007
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: apjaworski at mmm.com [mailto:apjaworski at mmm.com]
> Skickat: den 5 september 2007 18:24
> Till: Yu (Warren) Wang
> Kopia: r-help at stat.math.ethz.ch; r-help-bounces at stat.math.ethz.ch
> Ä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: apjaworski at mmm.com
> Tel: (651) 733-6092
> Fax: (651) 736-3122
>
>
>
> "Yu (Warren)
> Wang"
> <yu.wang at pdf.com> To
> Sent by: "r-help at stat.math.ethz.ch"
> r-help-bounces at st <r-help at stat.math.ethz.ch>
> 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
>
> ______________________________________________
> 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
> and provide commented, minimal, self-contained, reproducible code.
>
>
>
>
More information about the R-help
mailing list