#  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
  #
      print(llstuff$hs)  # print Hessian at end
  #
  # now get inverse of -Hessian for std errors
      backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) ) 
 #
 #  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
  #  get standard errors -- lazy, lazy
     solve(sol$hessian)
  #  done
    rm(list=ls())
  q()