[R] nls() and nls2() behavior?
ivo welch
ivowel at gmail.com
Tue May 11 16:46:26 CEST 2010
first, apologies for so many posts yesterday and today. I am
wrestling with nls() and nls2(). I have tried to whittle it down to a
simple example that still has my problem, yet can be cut-and-pasted
into R. here it is:
library(nls2)
options(digits=12);
y= c(0.4334,0.3200,0.5848,0.6214,0.3890,0.5233,0.4753,0.2104,0.3240,0.2827,0.3847,0.5571,0.5432,0.1326,0.3481)
x= c(0.3521,0.4334,0.3200,0.5848,0.6214,0.3890,0.5233,0.1379,0.2104,0.3240,0.3404,0.3847,0.5571,0.5432,0.1326)
dummy=c(rep( "m1", 7), rep( "m2", 3), rep( "m3", 5))
d= data.frame( y=y, x=x, dummy=dummy );
print(lm( y ~ as.factor(dummy) + x, data=d )) ## a simple fixed-effects model
objective.xs <- function( b, m1, m2, m3, x, ydebug ) {
tgts.mean= c(rep( m1, 7), rep( m2, 3), rep( m3, 5));
yhat= tgts.mean + b*x
#debug
cat("Parameters: b=", b, " m123=", m1, m2, m3, " leading to ",
sum( (ydebug-yhat)^2 ), "\n")
return( yhat );
}
sameformula= (y ~ objective.xs(b, m1, m2, m3, x, y))
cat("\nNLS2 Function --- what do you do?:\n");
print(nls2( sameformula, data=d, control= nls.control(tol=1e-12),
algorithm="brute-force",
start=list( b=0.3, m1=0.5, m2=0.5, m3=0.5, trace=TRUE ) ) );
cat("\n---but with this b starting value, this is so much better (?!):\n")
print(nls2( sameformula, data=d, control= nls.control(tol=1e-12),
algorithm="brute-force",
start=list( b=0.2, m1=0.5, m2=0.5, m3=0.5 ) ));
cat("\nNLS Function:\n");
print(nls( sameformula, data=d, control= nls.control(tol=1e-12),
start=list( b=0.3, m1=0.5, m2=0.5, m3=0.5, trace=TRUE ) ));
for some reason, nls2() does not start wandering off into a good
direction; and nls() has no direction (Error in nlsModel(formula, mf,
start, wts) : singular gradient matrix at initial parameter
estimates). I am probably doing something really silly, but I have
been staring at this example (not just the whittled down, but the
original one from which it was derived) for hours now, and I cannot
figure out how I have miscoded this. I hope the problem is obvious to
regular users...thanks for any helpful eyes here.
regards,
/iaw
----
Ivo Welch (ivo.welch at brown.edu, ivo.welch at gmail.com)
More information about the R-help
mailing list