[R] non-linear regression
Alane
citb332 at terra.com.br
Fri Dec 11 11:54:13 CET 2009
Katharine,
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.
The example below could be resolved in Excel, however in does not converge.
How to solve the problem?
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
07:36:00
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: tx.br.H.txt
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091211/840840e1/attachment-0002.txt>
More information about the R-help
mailing list