Keith Jewell
2009-Sep-21 16:17 UTC
[R] How to use nls when [selfStart] function returns NA or Inf??
Hi Everyone,
I posted this a couple of weeks ago with no responses. My interface (via
gmane) seemed a bit flakey at the time, so I'm venturing to repost with some
additional information.
I'm trying to write selfStart non-linear models for use with nls. In these
models some combinations of parameter values are illegal; the function value
is undefined.
That's OK when calling the function directly [e.g. SSmodel(x, pars...)]; it
is appropriate to return an appropriate non-value such as NA or Inf.
However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those
non-values lead to errors such as (but not limited to):
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
or (if I provide a gradient attribute)
Error in qr.default(.swts * attr(rhs, "gradient")) :
NA/NaN/Inf in foreign function call (arg 1)
A toy example demonstrating my problem (legal values of param are >1):
#-----------
SSexample<-selfStart(
model=function(x, param) x^log(param-1),
initial = function(mCall, data, LHS){
val<- 1.001
names(val) <- mCall[c("param")]
val
},
parameters=c("param")
)
#----------------
nls(y~SSexample(x, par), data=data.frame(x=1:10,y=rnorm(10)))
#---------
(repeat the last line a few times and you'll get the error).
I can't see a way of making nls either
a) stick to legal parameter values (which I'd have trouble specifying
anyway), or
b) (ideally) accept NA/NaN/Inf as indicating "bad" parameter values,
equivalent to very large errors in 'y' values
I really do want to use nls rather than a bounded optimisation tool (such as
optim) because this fits into a much bigger picture predicated on nls.
I'd appreciate any suggestions.
Keith Jewell
=================================sessionInfo() is given below.
I think the toy example above is enough to demonstrate my problem, but in
case it is relevant (I don't think it is) here is some more info about the
models I'm fitting:
I'm fitting sigmoidal models to microbial growth over time. The specific
model giving me problems at the moment is only one of a whole class of such
models which I need to work with. I specify it here to illustrate that it is
not always obvious what the bounds on the parameters are.
SSbaranyiBR94<-selfStart(
model=function(Time, y0, ymax, mu, lambda, m = 1, v = mu)
{
#+
# From the paper Baranyi J. & Roberts T. A. (1994). A dynamic approach to
predicting bacterial growth in food. Int J. of Fd. Micro. 23. 277-294
# Papers equations (6), (5b), (4b)
# eq 4b: y(Time) = y0 + mu' * A(Time) - ln(1+(exp(m' * mu' *
A(Time)) -
1)/exp(m' * (ymax-y0)))/m'
# eq 5b: A(Time) = Time + ln((exp(-v * Time) + q0)/(1+ q0))/v
# from eq 6: q0 = 1/(exp(mu'*lambda)-1)
#
# m' = m*ln(10)
# mu' = mu/ln(10)
#
# Cited paper gives defaults m = 1 and v = mu
#-
m <- m * log(10)
mu <- mu/log(10)
.value <- ymax - log(1 + (exp(m * (ymax - y0)) - 1)/exp(m * mu * (Time +
log((exp( - v * Time) + (1/(exp(v * lambda) - 1)))/(1 + (1/(exp(v *
lambda) - 1))))/v)))/m
.value
}
, initial = function(mCall, data, LHS)
{
##### <snip>
},
parameters=c("y0", "ymax", "mu",
"lambda")
)
===================> sessionInfo()
R version 2.9.1 (2009-06-26)
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.3 svSocket_0.9-43 svMisc_0.9-48 TinnR_1.0.3
R2HTML_1.59-1 Hmisc_3.6-1
loaded via a namespace (and not attached):
[1] cluster_1.12.0 grid_2.9.1 lattice_0.17-25 stats4_2.9.1
VGAM_0.7-9
Gabor Grothendieck
2009-Sep-21 16:23 UTC
[R] How to use nls when [selfStart] function returns NA or Inf??
With a small number of parameters just use brute force on grid to calculate starting values. See nls2 package. On Mon, Sep 21, 2009 at 12:17 PM, Keith Jewell <k.jewell at campden.co.uk> wrote:> Hi Everyone, > > I posted this a couple of weeks ago with no responses. My interface (via > gmane) seemed a bit flakey at the time, so I'm venturing to repost with some > additional information. > > I'm trying to write selfStart non-linear models for use with nls. In these > models some combinations of parameter values are illegal; the function value > is undefined. > That's OK when calling the function directly [e.g. ?SSmodel(x, pars...)]; it > is appropriate to return an appropriate non-value such as NA or Inf. > However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those > non-values lead to errors such as (but not limited to): > ? ?Error in numericDeriv(form[[3L]], names(ind), env) : > ? ? ?Missing value or an infinity produced when evaluating the model > or (if I provide a gradient attribute) > ? ?Error in qr.default(.swts * attr(rhs, "gradient")) : > ? ? ?NA/NaN/Inf in foreign function call (arg 1) > > A toy example demonstrating my problem (legal values of param are >1): > #----------- > SSexample<-selfStart( > ?model=function(x, param) x^log(param-1), > ?initial = function(mCall, data, LHS){ > ? val<- 1.001 > ? names(val) <- mCall[c("param")] > ? val > ? }, > ?parameters=c("param") > ) > #---------------- > nls(y~SSexample(x, par), data=data.frame(x=1:10,y=rnorm(10))) > #--------- > > (repeat the last line a few times and you'll get the error). > > I can't see a way of making nls either > a) stick to legal parameter values (which I'd have trouble specifying > anyway), or > b) (ideally) accept NA/NaN/Inf as indicating "bad" parameter values, > equivalent to very large errors in 'y' values > > I really do want to use nls rather than a bounded optimisation tool (such as > optim) because this fits into a much bigger picture predicated on nls. > > I'd appreciate any suggestions. > > Keith Jewell > =================================> sessionInfo() is given below. > > I think the toy example above is enough to demonstrate my problem, but in > case it is relevant (I don't think it is) here is some more info about the > models I'm fitting: > > I'm fitting sigmoidal models to microbial growth over time. The specific > model giving me problems at the moment is only one of a whole class of such > models which I need to work with. I specify it here to illustrate that it is > not always obvious what the bounds on the parameters are. > > SSbaranyiBR94<-selfStart( > model=function(Time, y0, ymax, mu, lambda, m = 1, v = mu) > { > ?#+ > ?# From the paper Baranyi J. & Roberts T. A. (1994). A dynamic approach to > predicting bacterial growth in food. Int J. of Fd. Micro. 23. 277-294 > ?# Papers equations (6), (5b), (4b) > ?# eq 4b: y(Time) = y0 + mu' * A(Time) - ln(1+(exp(m' * mu' * A(Time)) - > 1)/exp(m' * (ymax-y0)))/m' > ?# eq 5b: A(Time) = Time + ln((exp(-v * Time) + q0)/(1+ q0))/v > ?# from eq 6: q0 = 1/(exp(mu'*lambda)-1) > ?# > ?# ? m' = m*ln(10) > ?# ?mu' = mu/ln(10) > ?# > ?# Cited paper gives defaults m = 1 and v = mu > ?#- > ?m <- m * log(10) > ?mu <- mu/log(10) > ?.value <- ymax - log(1 + (exp(m * (ymax - y0)) - 1)/exp(m * mu * (Time + > ? ?log((exp( - v * Time) + (1/(exp(v * lambda) - 1)))/(1 + (1/(exp(v * > ? ?lambda) - 1))))/v)))/m > ? .value > } > , initial = function(mCall, data, LHS) > { > ##### <snip> > }, > parameters=c("y0", "ymax", "mu", "lambda") > ) > ===================>> sessionInfo() > R version 2.9.1 (2009-06-26) > 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.3 svSocket_0.9-43 ? ?svMisc_0.9-48 ? ? ?TinnR_1.0.3 > R2HTML_1.59-1 ? ? ?Hmisc_3.6-1 > > loaded via a namespace (and not attached): > [1] cluster_1.12.0 ?grid_2.9.1 ? ? ?lattice_0.17-25 stats4_2.9.1 > VGAM_0.7-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. >