[R] Fix for nls bug???
Keith Jewell
k.jewell at campden.co.uk
Tue Aug 5 18:19:05 CEST 2008
Hi All,
I've hit a problem using nls. I think it may be a restriction in the
applicability of nls and I may have found a fix, but I've been wrong before.
This example is simplified to the essentials. My real application is much
more complicated.
Take a function of matrix 'x' with additional arguments:
matrix 'aMat' whose values are _not_ to be determined by nls
vector 'Coeffs' whose vales _are_ to be determined.
For simplicity, this isn't a selfStart function with an 'initial' attribute,
but that doesn't change things.
Myfunc<-function(x, aMat, Coeffs)
{
#
# result = quadratic response in x with
# terms selected by aMat
#
aMat[aMat!=0] <- Coeffs
rowSums((x%*%aMat)%*%t(x))
}
If aMat is passed in by name (e.g. aMat = bMat) nls fails.
e.g.
#
# data frame with some noise
DF <- data.frame(x1 = runif(20, 1, 20), x2=runif(20, 1, 20))
DF$y <- 1 +DF$x1 +DF$x2 +DF$x1*DF$x2 +DF$x1^2 + DF$x2^2 + rnorm(20)
#
# matrix to pass in as aMat
bMat <- matrix(c(1,1,0,0), 2, 2)
#
# and nls fails
nls(y ~ Myfunc(cbind(x1, x2), bMat, aVec), DF, start=list(aVec=c(1,2)))
#
# pass in the same matrix other than by name and it works
nls(y ~ Myfunc(cbind(x1, x2), matrix(c(1,1,0,0), 2, 2), aVec), DF,
start=list(aVec=c(1,2)))
I think the problem lies in this line in nls
for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var), data,
env)
This adds values for some named arguments (bMat) as columns of the data
frame. The problem is that generally they don't have the same number of
rows. I've made it work for my example by replacing that line with this
line, which adds values for those arguments to the data frame as parameters
rather than as a column
attributes(mf)[["parameters"]] <-
c(attributes(mf)[["parameters"]],lapply(varNames[!varIndex], function(var)
eval(as.name(var), data, env)) )
Problem is, I really don't know nls internals enough to be sure I haven't
broken something.
And anyway, if this is really an improvement I ought to share it, but don't
know how.
Or I could have totally the wrong end of the stick...
Comments, corrections and advice are welcome.
Thanks in advance,
Keith Jewell
-----------------------
> sessionInfo()
R version 2.7.0 (2008-04-22)
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 svIO_0.9-5 R2HTML_1.58
[5] svMisc_0.9-5 svIDE_0.9-5
loaded via a namespace (and not attached):
[1] tools_2.7.0 VGAM_0.7-7
More information about the R-help
mailing list