Re: [R] question about non-linear least squares in R

From: Nilsson Fredrik X <Fredrik.X.Nilsson_at_skane.se>
Date: Tue, 11 Sep 2007 12:18:32 +0200


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.

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.