[R] question about nls

[Prof J C Nash (U30A)]

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")