[R] How to use nls when [selfStart] function returns NA or Inf??
Keith Jewell
k.jewell at campden.co.uk
Mon Sep 21 18:17:07 CEST 2009
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
More information about the R-help
mailing list