Numerical-MLE-for-Gamma-distsribution-Math-564.R

x<- c(13.9, 4.5, 7.7, 17.5, 7.4, 8.3, 20.9, 6.6, 1.6, 2, 27, 7.8, 3.5, 
      3.1, 7.6, 7.6, 5.1, 7.8, 11, 7.7)

loglike<- function(par, data){
  alpha <- par[1]
  beta <- par[2]
  n<- length(data)
  like<- -n*log(gamma(alpha)) - n* alpha * log(beta) + 
         (alpha - 1) * sum(log(data)) - sum(data)/beta
  return(like)
}

results<- optim(par = c(1, 1), data = x, fn = loglike, method = c("L-BFGS-B"),
                lower = rep(0.00000001, 2), control = list(fnscale = -1), 
                hessian = T)
results$par

## [1] 2.313558 3.859857

cov.est <- -1* solve(results$hessian)
cov.est

##            [,1]       [,2]
## [1,]  0.4698208 -0.7838323
## [2,] -0.7838323  1.6297004

results$par[1] - qnorm(1 - 0.05/2) * sqrt(cov.est[1, 1])

## [1] 0.970131

results$par[1] + qnorm(1 - 0.05/2) * sqrt(cov.est[1, 1])

## [1] 3.656986

results$par[2] - qnorm(1 - 0.05/2) * sqrt(cov.est[2, 2])

## [1] 1.357772

results$par[2] + qnorm(1 - 0.05/2) * sqrt(cov.est[2, 2])

## [1] 6.361941

score <- function(par, data){
  alpha <- par[1]
  beta <- par[2]
  n <- length(data)
  alpha.score <- -n * digamma(alpha) - n* log(beta) + sum(log(data))
  beta.score <- -alpha * n/beta + sum(data)/beta^2
  cbind(alpha.score, beta.score)
}

score(results$par, data = x)

##       alpha.score    beta.score
## [1,] 4.416895e-06 -5.184769e-06

info.fnt <- function(par, data){
  alpha <- par[1]
  beta <- par[2]
  n <- length(data)
  info.alpha <- -1 * n * trigamma(alpha)
  info.beta <- n * alpha/beta^2 - 2 * sum(data)/beta^3
  info.alpha.beta <- -1 * n/beta
  return(matrix(c(info.alpha, info.alpha.beta, info.alpha.beta,
                  info.beta = info.beta), nrow = 2))
}
info.fnt(results$par, data = x)

##            [,1]      [,2]
## [1,] -10.773166 -5.181539
## [2,]  -5.181539 -3.105759

results$hessian

##           [,1]     [,2]
## [1,] -10.77317 -5.18154
## [2,]  -5.18154 -3.10576

optim(par = c(1, 1), data = x, fn = loglike, method = c("L-BFGS-B"),
      lower = rep(0.00000001, 2), control = list(fnscale = -1), hessian = T,
      gr = score)

## $par
## [1] 2.313558 3.859857
## 
## $value
## [1] -60.57493
## 
## $counts
## function gradient 
##       15       15 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
## 
## $hessian
##           [,1]     [,2]
## [1,] -10.77317 -5.18154
## [2,]  -5.18154 -3.10576