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

Tue Sep 11 12:18:32 CEST 2007

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.