# augrain1.r                  
 #  Example 9.5 -- maximize likelihood/concentrated likelihood
 #
 #  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))
      unlike <-  n*theta[1]*log(theta[2]) + n*lgamma(theta[1]) - 
                    (theta[1]-1)*sumly + sumy/theta[2] }
 #
 #                  concentrated/profile log-likelihood function
    lcalph <- function (alpha) { # profile llike 
             n <- length(y) 
             sumy <- sum(y)
             sumly <- sum(log(y))
         lcalph <- - n*alpha*log(sumy/(n*alpha)) - n*lgamma(alpha) + 
                    (alpha-1)*sumly - n*alpha }
 #
 #                  define gradient function
    grunlike <- function (theta) { # gradient of unlike
             n <- length(y) 
             sumy <- sum(y)
             sumly <- sum(log(y))
             dldalph <- n*log(theta[2]) + n*digamma(theta[1]) - sumly
             dldbeta <- n*theta[1]/theta[2] - sumy/(theta[2]**2)
          grunlike <- c(dldalph,dldbeta) }
 #  
 #  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,v))
 #  maximize profile like
    that1 <- optimize(f=lcalph, c(a0/5,a0*5), maximum=TRUE )
    print("maximize concentrated loglikelihood")
    that1
    print("value at maximum")
    ta <-that1$maximum
    ta
    f <- lcalph(ta)
    f
 #  check gradient
    tb <- sy/(n*ta)
    print(c(ta,tb))
    grad <- grunlike(c(ta,tb))
    grad
 #  check original loglikelihood function
    f <- unlike( c(ta,tb) )
    f
 #                  done
    rm(list=ls())
  q()