# cinlls3.r                  November 2007, 2008, 2009
 #
 # simulation similar to Chapter 9, Exercise 24
 #
 # nonlinear least squares problem nllsqt1
 #  test problem from Fuller(1976)
 #
 # compare different methods for se's in nonlinear regression
 #  
  n <- 12
  p <- 3
  x <- c(0, .4, .8, 1.6, 0, 1, 2, 4, 0, 1.5, 3, 5.9) 
  x
 # no y's -- we'll be generating them
 #
 # function to compute everything: estimates, se's
 #
  methods <- function(y) {   # nonlinear least squares confidence int
 # create data frame or list for nls
  fulnls <- list(y,x)
 #  find better starting values
  t1 <- y[ x==max(x)]             # max value as asymptote
  xs <- x[ x < 2 ]                # use smallest x's
  ys <- y[ x < 2 ]
  ok <- ( t1 > ys )               # second subsetting
  xs <- xs[ok]
  ys <- log( t1 - ys[ok] )        # transform reversed y's
 #  simple linear regression
  t3 <- sum( (xs - mean(xs) )*ys )/ sum( (xs-mean(xs))**2 ) # slope
  t2 <- mean(ys)-t3*mean(xs)                                # intercept
 #  transform back from log and flip signs
  beta0 <- c(t1,-exp(t2),-t3)
 # sum of squares function
  Sbeta <- function(t) {      # sum of squares function
    e <- y - t[1] - t[2]*exp(-t[3]*x)
    Sbeta <- sum(e*e)  }
 # call nls
  soln <- nls( y ~ th1 + th2*exp(-th3*x),
        start=list(th1=beta0[1],th2=beta0[2],th3=beta0[3]) )
  beta <- coef(soln)
 # different methods for standard errors for coefficient estimates
    e <- y - beta[1] - beta[2]*exp(-beta[3]*x)    # residuals
  sse <- sum(e*e)                                 # error sum of squares
 # gradient matrix -- works as design matrix
  G <- matrix( c(rep(1,n),exp(-beta[3]*x),-beta[2]*x*exp(-beta[3]*x)),
          n,p)
 # meat for sandwich
  cs <- (n-1)*cov(G*e)
 # using other pieces -- first 2nd partial of gi * residuali
 #  start with zero and just do the nonzero terms
  del2ge <- matrix(0,p,p)
  del2ge[2,3] <- sum(e*(-x*exp(-beta[3])) )     # d2 g / db2 db3
  del2ge[3,2] <- del2ge[2,3]                    # symmetry
  del2ge[3,3] <- sum(e*(beta[2]*x*x*exp(-beta[3]*x))) # d2 g / db3 db3
  GG <- t(G) %*% G                      # analogue to X'X
  hdel2S <- GG - del2ge                 # half of hessian of S(beta)
 #
 # all of the different methods
  GGinv <- chol2inv(chol(GG))
  hdel2Sinv <- chol2inv(chol(hdel2S))
  m1 <- diag((sse/(n-p))*GGinv)           # sigma2hat*inv(G'G)
  m2 <- diag((sse/(n-p))*hdel2Sinv)       # using hessian
  m3 <- diag(GGinv%*%cs%*%GGinv)          # easy sandwich
  m4 <- diag(hdel2Sinv%*%cs%*%hdel2Sinv)  # harder sandwich
 # finish function by putting all together
  methods <- c(beta,sse,sqrt(c(m1,m2,m3,m4)))    } # end of function
 #
 # simulation experiment
  set.seed(5151917)                # set seed
  nrep <- 1000                     # number of replications
  betatrue <- c(150, -80, 1.5)     # true value of beta
  sigtrue <- 8                     # true noise level
 # mean vector is constant for all replications
  yhat <- betatrue[1] + betatrue[2]*exp(-betatrue[3]*x)
 #
 # generate Y's
  Ys <- matrix(rep(yhat,nrep)+sigtrue*rnorm(nrep*n),n,nrep)
 #
 # *** main step *** compute estimates, etc for each sample ***
  ests <- apply(Ys,2,methods)
 #
 # difference of estimates from true
  diffs <- ests[1:3,]-betatrue
 # sampling variation of estimates
  bias <- apply(diffs,1,mean)
  bias
  cov <- var(t(diffs))
 # is bias significant?
  sqrt(diag(cov)/nrep)        # statistically
  bias/betatrue               # practically
 # how do the four methods compare -- average se
  matrix(c(sqrt(diag(cov)),apply(ests[5:16,],1,mean)),3,5)
 #  sd(bhat) mean(se1) mean(se2) mean(se3) mean(se4)
 #
 #                  now look at coverage of confidence intervals
 #
  cover <- t(abs(matrix(t(diffs),nrep,12))) < 1.96*ests[5:16,]
  apply(cover,1,sum)
  matrix(apply(cover,1,sum),3,4)
  cover %*% t(cover)
 # compare method #1 with method #2 on first parameter
  comp <- rbind(abs(diffs[1,])<1.96*ests[5,],
                abs(diffs[1,])<1.96*ests[8,])
  apply(comp,1,sum)
  comp %*% t(comp)
 # done
  rm(list=ls())
  q()