I decided to follow up my own suggestion and look at the robustness of
nls vs. nlxb. NOTE: this problem is NOT one that nls() would usually be
applied to. The script below is very crude, but does illustrate that
nls() is unable to find a solution in >70% of tries where nlxb (a
Marquardt approach) succeeds.
I make no claim for elegance of this code -- very quick and dirty.
JN
debug<-FALSE
library(nlmrt)
x<-c(60,80,100,120)
y<-c(0.8,6.5,20.5,45.9)
mydata<-data.frame(x,y)
mydata
xmin<-c(0,0,0)
xmax<-c(8,8,8)
set.seed(123456)
nrep <- as.numeric(readline("Number of reps:"))
pnames<-c("a", "b", "d")
npar<-length(pnames)
# set up structure to record results
# need start, failnls, parnls, ssnls, failnlxb, parnlxb, ssnlxb
tmp<-matrix(NA, nrow=nrep, ncol=3*npar+4)
outcome<-as.data.frame(tmp)
rm(tmp)
colnames(outcome)<-c(paste("st-",pnames[[1]],''),
paste("st-",pnames[[2]],''),
paste("st-",pnames[[3]],''),
"failnls", paste("nls-",pnames[[1]],''),
paste("nls",pnames[[1]],''),
paste("nls-",pnames[[1]],''), "ssnls",
"failnlxb", paste("nlxb-",pnames[[1]],''),
paste("nlxb",pnames[[1]],''),
paste("nlxb-",pnames[[1]],''), "ssnlxb")
for (i in 1:nrep){
cat("Try ",i,":\n")
st<-runif(3)
names(st)<-pnames
if (debug) print(st)
rnls<-try(nls(y ~ exp(a + b*x)+d,start=st, data=mydata), silent=TRUE)
if (class(rnls) == "try-error") {
failnls<-TRUE
parnls<-rep(NA,length(st))
ssnls<-NA
} else {
failnls<-FALSE
ssnls<-deviance(rnls)
parnls<-coef(rnls)
}
names(parnls)<-pnames
if (debug) {
cat("nls():")
print(rnls)
}
rnlxb<-try(nlxb(y ~ exp(a + b*x)+d,start=st, data=mydata), silent=TRUE)
if (class(rnlxb) == "try-error") {
failnxlb<-TRUE
parnlxb<-rep(NA,length(st))
ssnlxb<-NA
} else {
failnlxb<-FALSE
ssnlxb<-rnlxb$ssquares
parnlxb<-rnlxb$coeffs
}
names(parnls)<-pnames
if (debug) {
cat("nlxb():")
print(rnlxb)
tmp<-readline()
cat("\n")
}
solrow<-c(st, failnls=failnls, parnls, ssnls=ssnls,
failnlxb=failnlxb, parnlxb, ssnlxb=ssnlxb)
outcome[i,]<-solrow
} # end loop
cat("Proportion of nls runs that failed = ",sum(outcome$failnls)/nrep,"\n")
cat("Proportion of nlxb runs that failed =
",sum(outcome$failnlxb)/nrep,"\n")
______________________________________________
R-help@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.
