>  # 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
 [1] 0.0 0.4 0.8 1.6 0.0 1.0 2.0 4.0 0.0 1.5 3.0 5.9
>  # 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
        th1         th2         th3 
 0.35920277 -0.45479742  0.05170761 
>   cov <- var(t(diffs))
>  # is bias significant?
>   sqrt(diag(cov)/nrep)        # statistically
       th1        th2        th3 
0.13459991 0.18562926 0.01100742 
>   bias/betatrue               # practically
        th1         th2         th3 
0.002394685 0.005684968 0.034471741 
>  # how do the four methods compare -- average se
>   matrix(c(sqrt(diag(cov)),apply(ests[5:16,],1,mean)),3,5)
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 4.2564228 4.2325788 4.2489003 3.4731201 3.5234068
[2,] 5.8701128 5.8755030 5.9127910 4.8915843 4.9787759
[3,] 0.3480853 0.3279204 0.3313207 0.2677572 0.2774447
>  #  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)
927 930 916 928 928 918 838 858 839 841 858 847 
>   matrix(apply(cover,1,sum),3,4)
     [,1] [,2] [,3] [,4]
[1,]  927  928  838  841
[2,]  930  928  858  858
[3,]  916  918  839  847
>   cover %*% t(cover)
 927 886 867 925 884 867 834 821 802 833 823 807
 886 930 861 887 926 863 807 857 790 810 855 797
 867 861 916 868 859 914 795 795 838 796 798 843
 925 887 868 928 885 868 836 822 805 836 824 810
 884 926 859 885 928 861 804 856 789 807 855 795
 867 863 914 868 861 918 795 797 839 796 799 846
 834 807 795 836 804 795 838 753 750 834 755 752
 821 857 795 822 856 797 753 858 730 756 852 737
 802 790 838 805 789 839 750 730 839 751 732 830
 833 810 796 836 807 796 834 756 751 841 758 753
 823 855 798 824 855 799 755 852 732 758 858 738
 807 797 843 810 795 846 752 737 830 753 738 847
>  # 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)
[1] 927 928
>   comp %*% t(comp)
     [,1] [,2]
[1,]  927  925
[2,]  925  928
>  # done
>   rm(list=ls())
>   q()