[R] non-linear regression

Katharine Mullen kate at few.vu.nl
Fri Dec 11 16:32:19 CET 2009 

>
> The problem of estimation of parameters in R is that you have to know the
> value of the initial estimates very accurately, otherwise it does not
> converge.

This was discussed on R-help in the last 2 weeks; see the thread on
'Starting estimates for nls Exponential Fit'.

>
> The example below could be resolved in Excel, however in  does not converge.
> How to solve the problem?

Can you consider a model with a different functional form or do you need
to fit this exact function?

If you want to minimize log(data) log(model) RSS, then:

tx.br <- read.table('tx.br.H.txt',header=F,dec=',')
tx.br <- tx.br[,1]
id <- 1:100

qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
  (A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))

newD <- log(tx.br)

HP <- nls(newD~log(qx.suav(id,A,B,C,D,E,F,G,H,K)),
          data=data.frame(id=id,newD=newD),
          trace=TRUE,
          nls.control(maxiter=50000,warnOnly=TRUE),
          algorithm='port',
          start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
            E=5.949082737,F=24.526811,G=0.000046733960,H=1.0970550987,
            K=0.771722501657),
          lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))

HP

par(mfrow=c(2,1))
matplot(cbind(fitted(HP), newD),type="l", main="model and fit")

matplot(cbind(exp(fitted(HP)), exp(newD)), type="l",
        main="transformed back to original space")

> I made the chart on a logarithmic scale to better visualize the differences.
>
> Send the data file attached.
>
> The commands are below:
>
> tx.br <- read.table('c:/tx.br.H.txt',header=F,dec=',')
> tx.br <-tx.br[,1]
> id<-1:100
>
> qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
>   (A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
>
> HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
>           data=data.frame(id=id,tx.br=tx.br),
>           trace=TRUE,nls.control(maxiter=50000,warnOnly=TRUE,minFactor =
> 0.1),
>          algorithm='port',
>          start=list(A=0.000644,B=0.016761290,C=0.10927095582,D=0.00094877,
>            E=5.949082737,F=24.526811,G=0.000046733960,H=1.0970550987,K=0.771722501657),
>           lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
>
> HP
>
> matplot(cbind(log(fitted(HP)), log(tx.br)),type="l")
>
>
>
> ----- Original Message -----
> From: "Katharine Mullen" <kate at few.vu.nl>
> To: "AneSR" <citb332 at terra.com.br>
> Cc: <r-help at r-project.org>
> Sent: Thursday, December 10, 2009 9:55 PM
> Subject: Re: [R] non-linear regression
>
>
> > You did not provide the data or a way of generating it.
> >
> > I would guess that Excel finds the same solution (the same residual sum-of
> > squares) as nls but that it uses a weak convergence criterion and/or does
> > not give you information regarding why it terminates.
> >
> > Regarding the step size:  you can set the minimum step size factor via the
> > minFactor argument of control.
> >
> > qx.suav <- function(id,A,B,C,D,E,F,G,H,K)
> >  (A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))
> >
> > ## make noisy data from model
> > id <- 1:1000
> > tx.br <- qx.suav(id,A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
> >                 E=0.0002127,F=38.5448420,G=0.0000115,H=1.1109286,
> >                 K=0.382070638)
> > set.seed(1)
> > tx.br <- tx.br + rnorm(length(tx.br),0,.0001)
> >
> > HP <- nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),
> >          data=data.frame(id=id,tx.br=tx.br),
> >          trace=TRUE,nls.control(maxiter=5000,warnOnly=TRUE),
> >          algorithm='port',
> >          start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,
> >            E=0.0002127,F=38.5448420,G=0.0000115,H=1.1109286,
> >            K=0.382070638),
> >          lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
> > matplot(cbind(fitted(HP), tx.br),type="l")
> >
> > On Thu, 10 Dec 2009, AneSR wrote:
> >
> >>
> >> I have a non-linear regression with 8 parameters to solve .... however it
> >> does not converge ... easily solves the excel ... including the initial
> >> estimates used in the R were found in the excel ... Another question is
> >> how
> >> to establish the increments of R by the parameters in the search ..
> >>
> >>
> >> qx.suav<-function(id,A,B,C,D,E,F,G,H,K){(A^((id+B)^C)+(D*exp(-E*(log(id)-log(F))^2))+(G*H^id)/(1+(K*G*H^id)))}
> >> HP<-nls(tx.br~qx.suav(id,A,B,C,D,E,F,G,H,K),data=data.frame(id=id,tx.br=tx.br),
> >> trace=TRUE,nls.control(maxiter=5000),algorithm='port',start=list(A=0.0006347,B=0.0453814,C=0.1353538,D=0.1353538,E=0.0002127,F=38.5448420,G=0.0000115,H=1.1109286,K=0.382070638),lower=list(A=0,B=0,C=0,D=0,E=0,F=0,G=0,H=0,K=0))
> >> summary(HP)
> >>
> >> How to solve this problem in R?
> >>
> >> Ane
> >> --
> >> View this message in context:
> >> http://n4.nabble.com/non-linear-regression-tp959977p959977.html
> >> Sent from the R help mailing list archive at Nabble.com.
> >>
> >> [[alternative HTML version deleted]]
> >>
> >> ______________________________________________
> >> 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.
> >>
> >
>
>
> --------------------------------------------------------------------------------
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
> Version: 8.5.426 / Virus Database: 270.14.102/2556 - Release Date: 12/10/09
> 07:36:00
>

More information about the R-help mailing list