# tests/fitGammaDist.R In EBukin/maxLik-dev: Maximum Likelihood Estimation and Related Tools

```## the idea and most commands were provided by Marco J. Maier, Institute for
## Statistics and Mathematics, Vienna University of Economics and Business

library(maxLik)
options(warn = -1, digits = 4 )
set.seed(5)
some_data <- rgamma(1e4, shape = 5, scale = 2)

# log-likelihood function(s)
logLL <- function(x, X)   # per observation for maxLik
dgamma(x = X, shape = exp(x[1]), scale = exp(x[2]), log = TRUE)
logLL_sum <- function(x, X)   # negative sum for nlm()
-sum(dgamma(x = X, shape = exp(x[1]), scale = exp(x[2]), log = TRUE))

sum(logLL(log(c(5,2)),some_data))
logLL_sum(log(c(5,2)),some_data)
all.equal( sum(logLL(log(c(5,2)),some_data)), -logLL_sum(log(c(5,2)),some_data))

d_logLL <- function(x, X){   # analytic 1. derivatives
cbind(shape=exp(x[1])*(-x[2]-psigamma(exp(x[1]),0)+log(X)),
scale= X / exp(x[2]) - exp(x[1]))
}

d_logLLNum <- function(x, X){
numericGradient( logLL, x, X = X )
}

colSums(d_logLL(log(c(5,2)),some_data))
colSums(d_logLLNum(log(c(5,2)),some_data))

all.equal( d_logLL(log(c(5,2)),some_data), d_logLLNum(log(c(5,2)),some_data),
check.attributes=FALSE)

# Hessian of log-likelihood function
dd_logLL <- function(x, X){   # analytic 2. derivatives
grad <- d_logLL( x, X )
hessian <- matrix(0, 2, 2)
hessian[1,1] <- sum( grad[,1] - exp(x[1])^2 * psigamma(exp(x[1]), 1) )
hessian[2,2] <- - sum( X / exp(x[2]) )
hessian[cbind(c(2,1), c(1,2))] <- -exp(x[1]) * length(X)
return(hessian)
}

dd_logLLNum <- function(x, X){
numericHessian( function(x,X) sum(logLL(x,X)), t0=x, X = X )
}
numericHessian( function(x,X) sum(logLL(x,X)),
grad = function(x,X) colSums(d_logLL(x,X)), x, X = X )
}

dd_logLL(log(c(5,2)),some_data)
dd_logLLNum(log(c(5,2)),some_data)
all.equal(dd_logLL(log(c(5,2)),some_data), dd_logLLNum(log(c(5,2)),some_data))
check.attributes=FALSE)

# estimation with nlm()
t_nlm <- system.time( r_nlm  <- nlm(logLL_sum, c(0,0), X=some_data, hessian=TRUE) )

# estimation with nlm() and gradients
result <- logLL_sum( x, X )
attr( result, "gradient" ) <- - colSums( d_logLL( x, X ) )
return( result )
}
t_nlmg <- system.time( r_nlmg  <- nlm(logLL_grad, c(0,0), X=some_data, hessian=TRUE) )

# estimation with nlm() and gradients and Hessian
logLL_hess <- function(x, X) {
result <- logLL_sum( x, X )
attr( result, "gradient" ) <- - colSums( d_logLL( x, X ) )
attr( result, "hessian" ) <- - dd_logLL( x, X )
return( result )
}
t_nlmgh <- system.time( r_nlmgh  <- nlm(logLL_hess, c(0,0), X=some_data, hessian=TRUE) )

# estimation with optim() / BFGS
t_bfgs <- system.time( r_bfgs <- optim(c(0,0), logLL_sum, X=some_data,
method="BFGS", hessian=TRUE) )

# estimation with maxLik() / BFGS
t_bfgsM <- system.time( r_bfgsM <- maxLik( logLL, start = c(0,0),
method="BFGS", X=some_data ) )

# estimation with maxLik() / BFGS with gradients
t_bfgsMg <- system.time( r_bfgsMg <- maxLik( logLL, d_logLL, start = c(0,0),
method="BFGS", X=some_data ) )

# estimation with maxLik() / BHHH
t_bhhh <- system.time( r_bhhh <- maxLik( logLL, start = c(0,0),
method="BHHH", X=some_data ) )

# estimation with maxLik() / BHHH with gradients
t_bhhhg <- system.time( r_bhhhg <- maxLik( logLL, d_logLL, start = c(0,0),
method="BHHH", X=some_data ) )

# estimation with maxLik() / NR
t_NRn <- system.time( r_NRn <- maxLik( logLL, start = c(0,0),
method="NR", X=some_data ) )

# estimation with maxLik() / NR with gradients
t_NRg <- system.time( r_NRg <- maxLik( logLL, d_logLL, start = c(0,0),
method="NR", X=some_data ) )

# estimation with maxLik() / NR with gradients and Hessian
t_NRgh <- system.time( r_NRgh <- maxLik( logLL, d_logLL, dd_logLL, start = c(0,0),
method="NR", X=some_data ) )

# log likelihood values
rbind(NLM=-r_nlm\$minimum,
BFGS=-r_bfgs\$value,
maxLikBfgs = logLik( r_bfgsM ),
BHHH = logLik( r_bhhh ),
NR_numeric= logLik( r_NRn ),

# estimated coefficients
pp <- exp(rbind(NLM=r_nlm\$estimate,
BFGS=r_bfgs\$par,
maxLikBfgs = coef( r_bfgsM ),
BHHH = coef( r_bhhh ),
NR_numeric= coef( r_NRn ),
colnames(pp) <- c("shape_alpha", "scale_theta")
pp

# some Hessians
-100*round(r_nlm\$hessian/100,0)
round(solve(r_nlm\$hessian),5)

-100*round(r_nlmg\$hessian/100,0)
round(solve(r_nlmg\$hessian),5)

-100*round(r_nlmgh\$hessian/100,0)
round(solve(r_nlmgh\$hessian),5)

-100*round(r_bfgs\$hessian/100,0)
round(solve(r_bfgs\$hessian),5)

100*round(r_NRn\$hessian/100,0)
round(solve(-r_NRn\$hessian),5)

100*round(r_NRg\$hessian/100,0)
round(solve(-r_NRg\$hessian),5)

# standard errors
se <- exp(rbind(NLM=sqrt(diag( solve(r_nlm\$hessian) )),
BFGS=sqrt(diag( solve(r_bfgs\$hessian) )),
maxLikBfgs = stdEr( r_bfgsM ),
BHHH = stdEr( r_bhhh ),