# augrain2.r                  
 #  Example 9.5 -- maximize likelihood/concentrated likelihood
 #    using nlm
 #
 #  data
  both <- matrix(scan("augrain.dat"),45,2,byrow=TRUE)
  y <- both[,2]
  y
 #                  define log-likelihood function
    unlike <- function (theta) { # loglikelihood for 2 param gamma
             n <- length(y) 
             sumy <- sum(y)
             sumly <- sum(log(y))
 #                  define gradient function
             dldalph <- n*log(theta[2]) + n*digamma(theta[1]) - sumly
             dldbeta <- n*theta[1]/theta[2] - sumy/(theta[2]**2)
      attr(unlike,"gradient") <- c(dldalph,dldbeta) 
 #                  define Hessian function
             d2ldaa <- n*trigamma(theta[1])
             d2ldab <- n/theta[2]
             d2ldbb <- -n*theta[1]/(theta[2]**2) + 2*sumy/(theta[2]**3)
      attr(unlike,"hessian") <- matrix(c(d2ldaa,d2ldab,d2ldab,d2ldbb),2,2) 
 #                  function
      unlike <-  n*theta[1]*log(theta[2]) + n*lgamma(theta[1]) - 
                    (theta[1]-1)*sumly + sumy/theta[2] }
 #  
 #                  define log-likelihood function without grad, hess
    unlike0 <- function (theta) { # loglikelihood for 2 param gamma
             n <- length(y) 
             sumy <- sum(y)
             sumly <- sum(log(y))
 #                  function
      unlike0 <-  n*theta[1]*log(theta[2]) + n*lgamma(theta[1]) - 
                    (theta[1]-1)*sumly + sumy/theta[2] }
 #  
 #  starting values from method of moments
    n <- length(y)
    sy <- sum(y)
    vy <- var(y)
    b0 <- n*vy/sy
    a0 <- sy*sy/(n*n*vy)
    v <- unlike(c(a0,b0))
    print("method of moments estimates for starting values")
    print(c(a0,b0))
    print(v)
 #
 #                               now use nlm with gradient, Hessian
    that1 <- nlm(unlike, c(a0,b0), hessian=TRUE, typsize=c(1,1))
    that1
 #
 #                               now use nlm without gradient, Hessian
    that1 <- nlm(unlike0, c(a0,b0), hessian=TRUE, 
                 print.level=1, typsize=c(1,1))
    that1
 #
 #  get inverse of Hessian for std errors
    solve(that1$hessian)        # or should I use Cholesky?
 #
 #                  done
    rm(list=ls())
  q()