>  # oranget3.r                  October 2008
>  #
>  # fit orange tree data -- just third tree
>  #  from Draper & Smith
>  # 
>   x <- c(118,484,664,1004,1231,1372,1582)
>   y <- c(30,51,75,108,115,139,140)
>  #
>  # simple linear regression
>   slr <- function(y,x) {
+     xb <- mean(x)
+     yb <- mean(y)
+     b1 <- sum((x-xb)*(y-yb))/sum((x-xb)*(x-xb))
+     b0 <- yb - b1*xb
+     slr <- c(b0,b1)}
>  #
>  # find starting values
>  #
>   ymx <- max(y)
>  #
>  # get other starting values
>   Inot <- (y != ymx)
>   Inot
[1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
>   b0b1 <- slr( log(y[Inot]/(ymx-y[Inot])),x[Inot])
>   b0b1
[1] -2.361694845  0.004131767
>   third <- list(y,x)
>  # 
>  # call nonlinear least squares code nls
>  #
>   sol1 <- nls(y~th1/(1+exp(th2-th3*x)),data=third,
+    start=list(th1=ymx,th2=-b0b1[1],th3=b0b1[2]),trace=T )
668.8307 :  1.400000e+02 2.361695e+00 4.131767e-03 
668.1247 :  1.462007e+02 1.596177e+00 2.267760e-03 
156.8811 :  1.593868e+02 1.841021e+00 2.486950e-03 
156.6709 :  1.588209e+02 1.832690e+00 2.494111e-03 
156.6708 :  1.588295e+02 1.832730e+00 2.494046e-03 
>   sol1
Nonlinear regression model
  model:  y ~ th1/(1 + exp(th2 - th3 * x)) 
   data:  third 
      th1       th2       th3 
1.588e+02 1.833e+00 2.494e-03 
 residual sum-of-squares: 156.7

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.173e-06 
>  # covariance matrix of coefficient estimates
>   vcov(sol1)
              th1           th2           th3
th1 227.561826414 -9.195143e-01 -6.230956e-03
th2  -0.919514312  3.849845e-02  6.019525e-05
th3  -0.006230956  6.019525e-05  2.137668e-07
>  # 
>  # second way using simple Gauss-Newton
>  # 
>   tree <- function(theta) {       # find g, G, and e
+      g <- theta[1]/(1 + exp(theta[2]-theta[3]*x))
+      e <- y - g
+      G <- matrix(0,7,3)
+      G[,1] <- 1/(1+exp(theta[2]-theta[3]*x))
+      G[,2] <- -g*exp(theta[2]-theta[3]*x)/(1+exp(theta[2]-theta[3]*x))
+      G[,3] <- g*exp(theta[2]-theta[3]*x)*x/(1+exp(theta[2]-theta[3]*x)) 
+      tree <- list(e,g,G) }
>  #
>    beta <- c(ymx,-b0b1[1],b0b1[2])       # starting values
>    beta
[1] 1.400000e+02 2.361695e+00 4.131767e-03
>    for ( k in(1:8) ) {                  # loop start on k
+      treek <- tree(beta)                # get e,g,G
+      sse <- sum(treek[[1]]^2)           # sse
+      print("k,SSE")
+      print(c(k,sse))
+      G <- treek[[3]]
+      GG <- t(G) %*% G
+      L <- t( chol(GG) )                 # factor of G'G
+      step <- backsolve(t(L),forwardsolve(L, t(G) %*% treek[[1]] ) )
+      beta <- beta + step                # take step
+      print("beta,step")
+      print(c(beta,step))      }         # end of loop on k
[1] "k,SSE"
[1]   1.0000 668.8307
[1] "beta,step"
[1] 146.200704542   1.596177355   0.002267760   6.200704542  -0.765517490
[6]  -0.001864007
[1] "k,SSE"
[1]   2.0000 668.1245
[1] "beta,step"
[1] 1.593868e+02 1.841021e+00 2.486950e-03 1.318610e+01 2.448434e-01
[6] 2.191891e-04
[1] "k,SSE"
[1]   3.0000 156.8811
[1] "beta,step"
[1]  1.588209e+02  1.832690e+00  2.494111e-03 -5.658857e-01 -8.331058e-03
[6]  7.161809e-06
[1] "k,SSE"
[1]   4.0000 156.6709
[1] "beta,step"
[1]  1.588295e+02  1.832730e+00  2.494046e-03  8.557753e-03  4.073703e-05
[6] -6.573397e-08
[1] "k,SSE"
[1]   5.0000 156.6708
[1] "beta,step"
[1]  1.588295e+02  1.832730e+00  2.494046e-03 -2.210149e-05 -2.311933e-07
[6]  3.365500e-10
[1] "k,SSE"
[1]   6.0000 156.6708
[1] "beta,step"
[1]  1.588295e+02  1.832730e+00  2.494046e-03  1.353944e-07  1.048402e-09
[6] -2.365313e-12
[1] "k,SSE"
[1]   7.0000 156.6708
[1] "beta,step"
[1]  1.588295e+02  1.832730e+00  2.494046e-03 -4.057834e-10 -1.091269e-11
[6] -4.277951e-16
[1] "k,SSE"
[1]   8.0000 156.6708
[1] "beta,step"
[1]  1.588295e+02  1.832730e+00  2.494046e-03  1.050873e-11 -7.071790e-14
[6] -3.308997e-16
>  #  get covariance matrix -- divide by n or n-p ???
>      Covb <- (sse/4)*backsolve(t(L),forwardsolve(L,diag(c(1,1,1))) )
>      Covb
              [,1]          [,2]          [,3]
[1,] 227.561591966 -9.195146e-01 -6.230954e-03
[2,]  -0.919514644  3.849846e-02  6.019529e-05
[3,]  -0.006230954  6.019529e-05  2.137669e-07
>  # done
>   rm(list=ls())
>  q()