From: Nilsson Fredrik X <Fredrik.X.Nilsson_at_skane.se>

Date: Tue, 11 Sep 2007 12:18:32 +0200

# I do this a couple of times until convergence

apa3<-optim(apa3$par, foqufu, control=list(maxit=20000, reltol=1e-10))

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

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.

R-help_at_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. Received on Tue 11 Sep 2007 - 10:27:23 GMT

Date: Tue, 11 Sep 2007 12:18:32 +0200

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)

}

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:

- Scale your y's by dividing all values by 1e6.
- Plot x vs. y. The plot looks like a quadratic function.
- Fit a quadratic const. + B*x^2 - this a linear regression problem so use lm.
- Plot the predictions.
- Eyball the necessary shift - MA is around 0.01. Refit const. + B*(x-.01)^2. Should get const.=1.147 and B=139.144
- 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.

R-help_at_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. Received on Tue 11 Sep 2007 - 10:27:23 GMT

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

Archive generated by hypermail 2.2.0, at Wed 12 Sep 2007 - 02:35:38 GMT.

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