[R] nls with plinear and function on RHS
Keith Jewell
k.jewell at campden.co.uk
Thu Oct 2 17:54:01 CEST 2008
Wonderful, great! That solves my problem nicely (on to the next one!).
Thanks a lot,
Keith Jewell
"Katharine Mullen" <kate at few.vu.nl> wrote in message
news:Pine.GSO.4.56.0810021616340.138 at laurel.few.vu.nl...
> a) you should put values for Ca, Cb, Cc directly into the data list as
> data=list(Ca=1, ....
>
> b) you can simplify the call to
>
> # idealised data set
> aDF <- data.frame( x= c(1.80, 9.27, 6.48, 2.61, 9.86, 5.93, 6.76, 5.52,
> 6.06, 8.62),
> y= c(24.77, 2775.07, 895.15, 60.73, 3373.57, 677.82, 1021.92, 542.84,
> 725.25, 2200.04))
>
> bFunc <- function(x, Cd) cbind(Ca=1,Cb=x, Cc=x^Cd)
>
> # nls, plinear algorithm, RHS from function
> nls(y ~ bFunc(x, Cd), data=list(x=aDF$x, y=aDF$y),
> start=list(Cd=3), algorithm="plinear")
>
> On Thu, 2 Oct 2008, Keith Jewell wrote:
>
>> Dear R gurus,
>>
>> As part of finding initial values for a much more complicated fit I want
>> to
>> fit a function of the form y ~ a + bx + cx^d to fairly "noisy" data and
>> have
>> hit some problems.
>>
>> To demonstrate the specific R-related problem, here is an idealised data
>> set, smaller and better fitting than reality:
>> # idealised data set
>> aDF <- data.frame( x= c(1.80, 9.27, 6.48, 2.61, 9.86, 5.93, 6.76, 5.52,
>> 6.06, 8.62),
>> y= c(24.77, 2775.07, 895.15, 60.73, 3373.57, 677.82, 1021.92, 542.84,
>> 725.25, 2200.04))
>>
>> And here are some starting values, far better than I'd have in reality
>> # good starting values
>> startL <- list(Ca=4, Cb=3, Cc=3, Cd=3)
>>
>> In these idealised circumstances nls converges using the default
>> algorithm.
>> # nls, default algorithm
>> nls(y ~ Ca + Cb*x + Cc*x^Cd, data=aDF, start=startL)
>>
>> Unfortunately, in reality it often fails to converge. This model is
>> linear
>> in a, b and c so I've used the plinear algorithm.
>> # nls, plinear algorithm, explicit RHS
>> nls(y ~ cbind(Ca=1,Cb=x, Cc=x^Cd), data=aDF, start=startL["Cd"],
>> algorithm="plinear")
>>
>> This converges much more often but sometimes crashes with the error
>> message
>> "Error in numericDeriv(form[[3]], names(ind), env) :
>> Missing value or an infinity produced when evaluating the model"
>>
>> I deduce that it is failing in the numerical differentiation of x^Cd (but
>> don't know why), so I thought I'd avoid the numerical differentiation by
>> putting the RHS in a function to which I could (later) add a 'gradient'
>> attribute
>> # function to provide RHS
>> bFunc <- function(x, Ca, Cb, Cc, Cd) cbind(Ca=1,Cb=x, Cc=x^Cd)
>>
>> # nls, plinear algorithm, RHS from function
>> nls(y ~ bFunc(x, Ca, Cb, Cc, Cd), data=aDF, start=startL["Cd"],
>> algorithm="plinear")
>>
>> However, this gives me
>> "Error in nls(y ~ bFunc(x, Ca, Cb, Cc, Cd), data = aDF, start =
>> startL["Cd"], :
>> parameters without starting value in 'data': Ca, Cb, Cc"
>>
>> Can anyone tell me
>> a) why putting the RHS into a function "broke" the plinear algorithm
>> b) if there's a better approach to my problem
>>
>> Thanks in advance,
>>
>> Keith Jewell
>> -----------------
>> I'm using V2.7.2...
>> > sessionInfo()
>> R version 2.7.2 (2008-08-25)
>> i386-pc-mingw32
>>
>> locale:
>> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
>> Kingdom.1252;LC_MONETARY=English_United
>> Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>>
>> attached base packages:
>> [1] stats graphics grDevices datasets tcltk utils methods
>> base
>>
>> other attached packages:
>> [1] xlsReadWrite_1.3.2 svSocket_0.9-5 TinnR_1.0.2 R2HTML_1.59
>> Hmisc_3.4-3
>>
>> loaded via a namespace (and not attached):
>> [1] cluster_1.11.11 grid_2.7.2 lattice_0.17-14 svMisc_0.9-5
>> VGAM_0.7-7
>>
>> ... but have also tried todays V2.7.2 patched and V2.8.0alpha, both of
>> which
>> give the same behaviour
>>
>> ______________________________________________
>> 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.
>>
>
> ______________________________________________
> 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.
>
More information about the R-help
mailing list