>  #  chex94.r        ML analysis of Extended Example
>  #                  October 2007, 2008
>  #
>  # first define likelihood function
>     lglik <- function(t) { # log-likelihood for extended example
+         t1 <- t[1]         #  includes gradient and Hessian
+         t2 <- t[2]
+         n <- c(17,182,60,176)           # cell counts
+         p <- c( 2*t1*t2, t1*(2-t1-2*t2), t2*(2-t2-2*t1), (1-t1-t2)**2 )
+         ll <- crossprod( n, log(p) )    # likelihood
+           # Jacobian matrix of probs (4) by params (2)
+         del <- matrix( c(2*t2,2*t1,2*(1-t1-t2),-2*t1,
+           -2*t2,2*(1-t1-t2),-2*(1-t1-t2),-2*(1-t1-t2) ),2,4)
+         gradient <- del %*% (n/p)       # gradient vector
+         del2 <- array( c(0,2,2,0, -2,-2,-2,0, 0,-2,-2,-2, 2,2,2,2),
+                          c(2,2,4) ) 
+         hessian <-            # Hessian matrix
+            n[1]*( del2[,,1]/p[1] - outer(del[,1]/p[1],del[,1]/p[1]) ) +
+            n[2]*( del2[,,2]/p[2] - outer(del[,2]/p[2],del[,2]/p[2]) ) +
+            n[3]*( del2[,,3]/p[3] - outer(del[,3]/p[3],del[,3]/p[3]) ) +
+            n[4]*( del2[,,4]/p[4] - outer(del[,4]/p[4],del[,4]/p[4]) ) 
+           # return a list of information
+           lglik <- list(f=ll,gr=gradient,hs=hessian)  }
>   ##      lglik 
>   #
>   #  hope this starting point works
>     t <- c(.263,.074)
>   #
>     for (i in 1:5) {     # loop for Newton iterations
+       llstuff <- lglik(t)
+       print("iteration")
+       print(i)
+       print("function, gradient")
+       print(c(llstuff$f,llstuff$gr))
+       L <- t(chol(-llstuff$hs))         # note negative
+       step <- backsolve( t(L),forwardsolve(L,llstuff$gr) )
+       t <- t + step      # already negated Hessian
+       print("current params, step")
+       print(c(t,step))  
+                    }     # end iteration loop
[1] "iteration"
[1] 1
[1] "function, gradient"
[1] -494.6769   25.4822  237.6884
[1] "current params, step"
[1] 0.265556386 0.089487453 0.002556386 0.015487453
[1] "iteration"
[1] 2
[1] "function, gradient"
[1] -492.6031294   -0.4043816   37.1720467
[1] "current params, step"
[1]  0.264479789  0.093036913 -0.001076597  0.003549460
[1] "iteration"
[1] 3
[1] "function, gradient"
[1] -4.925354e+02  2.462549e-03  1.290454e+00
[1] "current params, step"
[1]  2.644444e-01  9.316864e-02 -3.542888e-05  1.317306e-04
[1] "iteration"
[1] 4
[1] "function, gradient"
[1] -4.925353e+02 -1.521763e-06  1.639368e-03
[1] "current params, step"
[1]  2.644443e-01  9.316881e-02 -4.635360e-08  1.679047e-07
[1] "iteration"
[1] 5
[1] "function, gradient"
[1] -4.925353e+02 -6.746825e-13  2.656073e-09
[1] "current params, step"
[1]  2.644443e-01  9.316881e-02 -7.462867e-14  2.719867e-13
>   #
>       print(llstuff$hs)  # print Hessian at end
          [,1]       [,2]
[1,] -3900.903  -1067.863
[2,] -1067.863 -10058.455
>   #
>   # now get inverse of -Hessian for std errors
>       backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) ) 
              [,1]          [,2]
[1,]  2.640242e-04 -2.803030e-05
[2,] -2.803030e-05  1.023947e-04
>  #
>  #  do again with nlm
>  #
>  #  define likelihood function
>     ull <- function(t) { # un-log-likelihood for extended example
+         t1 <- t[1]         #  without gradient and Hessian
+         t2 <- t[2]
+         n <- c(17,182,60,176)             # cell counts
+         p <- c( 2*t1*t2, t1*(2-t1-2*t2), t2*(2-t2-2*t1), (1-t1-t2)**2 )
+         ull <- -crossprod( n, log(p) ) }  # unlikelihood
>   #  minimize using nlm
>      sol <- nlm(ull, c(.263,.074), hessian=T )
>      sol
$minimum
[1] 492.5353

$estimate
[1] 0.26444429 0.09316885

$gradient
[1] -6.392298e-05  4.953620e-05

$hessian
         [,1]      [,2]
[1,] 3899.064  1068.172
[2,] 1068.172 10039.791

$code
[1] 1

$iterations
[1] 8

>   #  get standard errors -- lazy, lazy
>      solve(sol$hessian)
              [,1]          [,2]
[1,]  2.641716e-04 -2.810623e-05
[2,] -2.810623e-05  1.025940e-04
>   #  done
>     rm(list=ls())
>   q()