[R] How can I fitting this function from values.

[Rolf Turner]

On 13/11/2007, at 8:54 PM, Frede Aakmann Tøgersen wrote:

> R is case sensitive ;-) so in your list of start values you must  
> have B=3 and not b=3.

	Duhhhhhh!!! Typo.

> Perhaps you have an object named b of a different length than your  
> data somewhere giving rise to your error message.

	Nope.  That's not it.
>
> Also brute force wont work as Joerg point out. We have
>
> A*1/(1+t/B)*1/sqrt(1+t/C) = D/(B+t)/sqrt(C+t),
>
> where D = A*B*sqrt(C).
> We notice that we need B > -t and C > -t.

OTOH if one invokes optim() --- with method="BFGS" (no constraints)
it goes like a train (given the excellent starting values, at least):

foo <- function(t,ccc){
ccc[1]*1/(1+t/ccc[2])*1/sqrt(1+1/ccc[3])
}
ttt <- seq(1,10,length=100)
yyy <- foo(ttt,c(5,3,2))
set.seed(42)
dat.test <- data.frame(t=ttt,y=yyy+rnorm(100,0,0.1))
with(dat.test,plot(t,y))
with(dat.test,lines(t,foo(t,c(5,3,2))))

bar <- function(par,y,t) {
         yhat <- foo(t,par)
         sum((y-yhat)^2)
}

fit.test <- optim(c(5,3,2),bar,method="BFGS",control=list 
(trace=6),t=dat.test$t,y=dat.test$y)

It converges virtually instantaneously and the fit looks good:

 > fit.test
$par
[1] 5.096597 2.888116 2.039852

$value
[1] 1.062285

$counts
function gradient
       17        6

$convergence
[1] 0

$message
NULL

 > yhat <- foo(dat.test$t,fit.test$par)
 > lines(dat.test$t,yhat,col="red")

Moral of the story:  It's better and easier and less frustrating to  
use optim()
rather than nls() wherever possible.

		cheers,

			Rolf Turner
>
> Now one can use the "port" algorithme of nls() with constraints on  
> parameters:
>
>
> foo <- function(t,ccc){
>   D <- ccc[1]
>   B <- ccc[2]
>   C <- ccc[3]
>   D/(B+t)/sqrt(C+t)
> }
>
> ttt <- seq(1,10,length=100)
> yyy <- foo(ttt,c(5*3*sqrt(2),3,2))
> set.seed(42)
>
> dat.test <- data.frame(t=ttt,y=yyy+rnorm(100,0,.1))
>
>
> fit.test <- nls(y~D/(B+t)/sqrt(C+t),data=dat.test,
>                 start=c(D=5*3*sqrt(2),B=3,C=2),
> 		    lower=c(-Inf,-1,-1),alg="port")
> # -Inf could be 0, -1 could be -min(dat.test$t)
>
> with(dat.test,plot(t,y))
> with(dat.test,lines(t,foo(t,c(21.21,3,2)),col="green")) # Looks cool.
> lines(ttt,fitted(fit.test),col="red")
>
>
> Best regards
>
> Frede Aakmann Tøgersen
> Scientist
>
>
> UNIVERSITY OF AARHUS
> Faculty of Agricultural Sciences
> Dept. of Genetics and Biotechnology
> Blichers Allé 20, P.O. BOX 50
> DK-8830 Tjele
>
> Phone:   +45 8999 1900
> Direct:  +45 8999 1878
>
> E-mail:  FredeA.Togersen at agrsci.dk
> Web:	   http://www.agrsci.org				
>
> This email may contain information that is confidential.
> Any use or publication of this email without written permission  
> from Faculty of Agricultural Sciences is not allowed.
> If you are not the intended recipient, please notify Faculty of  
> Agricultural Sciences immediately and delete this email.
>
>
>
>
>> -----Oprindelig meddelelse-----
>> Fra: r-help-bounces at r-project.org
>> [mailto:r-help-bounces at r-project.org] På vegne af Rolf Turner
>> Sendt: 12. november 2007 21:16
>> Til: Horacio Castellini
>> Cc: r-help at r-project.org
>> Emne: Re: [R] How can I fitting this function from values.
>>
>>
>> On 13/11/2007, at 8:27 AM, Horacio Castellini wrote:
>>
>>> Hi all, I'd like fit this function:
>>>
>>> G(t)=A*1/(1+t/B)*1/sqrt(1+t/C)
>>>
>>> where A, B and C are fitting constants that one like search. I have
>>> got a fcs<-(t,G(t)) value table in a data frame.
>>>
>>> Any one can help me? Tahnks Horacio.
>>
>>
>> I ***thought*** that nls would solve your problem in a trice.
>> But I tried a toy example before replying to you, and it messed up
>> mightily:
>>
>> foo <- function(t,ccc){
>> ccc[1]*1/(1+t/ccc[2])*1/sqrt(1+1/ccc[3])
>> }
>> ttt <- seq(1,10,length=100)
>> yyy <- foo(ttt,c(5,3,2))
>> set.seed(42)
>> dat.test <- data.frame(t=ttt,y=yyy+rnorm(100,0,0.1))
>> with(dat.test,plot(t,y))
>> with(dat.test,lines(t,foo(t,c(5,3,2)))) # Looks cool.
>> fit.test <- nls(y~A*1/(1+t/B)*1/sqrt(1+t/C),start=list
>> (A=5,b=3,C=2),data=dat.test)
>> # Gives error; no idea what the implications are.
>>
>> The error message:
>>
>> Error in `[[<-.data.frame`(`*tmp*`, var, value =
>> c(0.70934153524875, 0.761288405463172,  :
>>    replacement has 65 rows, data has 100 In addition: Warning  
>> message:
>> In t/B : longer object length is not a multiple of shorter
>> object length
>>
>> I don't think I have *ever* had a good experience with nls();
>> it is always pretty flaky.
>> (Well, maybe it worked for me once ..., can't remember for  
>> sure. :-) )
>>
>> 		cheers,
>>
>> 			Rolf Turner
>>
>>
>> ##################################################################### 
>> #
>> Attention:\ This e-mail message is privileged and
>> confid...{{dropped:9}}
>>
>> ______________________________________________
>> 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.
>>

######################################################################
Attention: 
This e-mail message is privileged and confidential. If you are not the 
intended recipient please delete the message and notify the sender. 
Any views or opinions presented are solely those of the author.

This e-mail has been scanned and cleared by MailMarshal 
www.marshalsoftware.com
######################################################################