R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  # 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
 [1] 4.63 4.18 6.69 1.77 3.66 1.37 3.65 1.95 1.88 5.86 6.45 2.53 5.07 2.94 1.37
[16] 4.05 0.92 2.20 4.95 4.15 3.82 2.85 3.36 4.68 5.34 3.99 7.67 3.84 4.74 3.94
[31] 4.20 4.71 2.33 4.28 4.20 5.89 5.43 3.95 2.98 1.19 2.85 3.88 6.97 2.89 2.60
>  #                  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")
[1] "method of moments estimates for starting values"
>     print(c(a0,b0))
[1] 5.8357604 0.6582023
>     print(v)
[1] 84.7093
>  #
>  #                               now use nlm with gradient, Hessian
>     that1 <- nlm(unlike, c(a0,b0), hessian=TRUE, typsize=c(1,1))
>     that1
$minimum
[1] 84.50925

$estimate
[1] 5.1414032 0.7470934

$gradient
[1] 2.321763e-07 6.536993e-07

$hessian
          [,1]     [,2]
[1,]  9.657366  60.2294
[2,] 60.229399 414.2976

$code
[1] 1

$iterations
[1] 7

>  #
>  #                               now use nlm without gradient, Hessian
>     that1 <- nlm(unlike0, c(a0,b0), hessian=TRUE, 
+                  print.level=1, typsize=c(1,1))
iteration = 0
Step:
[1] 0 0
Parameter:
[1] 5.8357604 0.6582023
Function Value
[1] 84.7093
Gradient:
[1] 0.5522702675 0.0003030323

iteration = 7
Parameter:
[1] 5.1414032 0.7470934
Function Value
[1] 84.50925
Gradient:
[1] 2.321763e-07 6.536993e-07

Relative gradient close to zero.
Current iterate is probably solution.

>     that1
$minimum
[1] 84.50925

$estimate
[1] 5.1414032 0.7470934

$gradient
[1] 2.321763e-07 6.536993e-07

$hessian
          [,1]     [,2]
[1,]  9.657366  60.2294
[2,] 60.229399 414.2976

$code
[1] 1

$iterations
[1] 7

>  #
>  #  get inverse of Hessian for std errors
>     solve(that1$hessian)        # or should I use Cholesky?
           [,1]        [,2]
[1,]  1.1094002 -0.16128143
[2,] -0.1612814  0.02586036
>  #
>  #                  done
>     rm(list=ls())
>   q()