## READ DATA from file!

sessionInfo()


OrangeTreeData.y <- t(read.table("../DATA/OrangeTree.dat",skip=7))
OrangeTreeData.x<-scan(file='../DATA/OrangeTree.dat',what=numeric(), nlines=1, skip=5)

treedata<-data.frame(Tree=as.vector(col(OrangeTreeData.y)), age=rep(OrangeTreeData.x,5), circumference=as.vector(OrangeTreeData.y))

#t1 <- system.time(nfit <- try(nls(circumference ~ phi1 / (1 + exp(-(age-phi2)/phi3)), data=Orange, start = list(phi1 = 200, phi2=800, phi3=400))))["elapsed"]
t1 <- system.time(nfit <- try(nls(circumference ~ phi1 / (1 + exp(-(age-phi2)/phi3)), data=treedata,
                                  start = list(phi1 = 200, phi2=800, phi3=400))))["elapsed"]


if (class(nfit)=="try-error")  {  # or if (inherits(nfit,"try-error"))
  results <- list(obj=NA,coef=rep(NA,4),
                  sd=rep(NA,4))
} else {
  est.sd <- sqrt(deviance(nfit)/(nobs(nfit)-length(coef(nfit))))
  logSigma <- log(est.sd)
  ## can we get uncertainty for logSigma??
  results <- list(obj=-logLik(nfit),  ## objective function/Neg log likelihood (NULL for MCMC)
                coef=c(coef(nfit),logSigma=logSigma),  ## coefficients
                  sd=sqrt(diag(vcov(nfit))),
                  vcov=vcov(nfit),
                  confint.quad=confint.default(nfit),
                  confint.profile=confint(nfit),
                  time=c(fit=unname(t1))
                  )
}

save("results",file="fit.RData")