Parameter estimation of Gamma Distribution using R

Question

This question might also be suited to the programming site, but I thought since there is enough on the statistics side, I would use this forum. I am trying to estimate the alpha parameter in a Gamma distribution using maximum likelihood method, and using the optimization functions available in R.

To begin with, I generated a random sample from Gamma(Alpha, Beta) in R.

shape <- 2
scale <- 1.5
set.seed(123456)
myData <- round(rgamma(n=50, shape=shape, scale=scale),2)

Using the maximum likelihood estimation method, and setting up the likelihood function to be in terms of alpha only, I created a function in R and I am trying to optimize it. So I wrote the likelihood function, took the log, took the partial derivative with respect to Beta, and found the MLE of Beta. I then substituted the MLE of Beta back into the likelihood function to arrive at the likelihood in terms of alpha only. My function is as follows:

objFunction <- function(myData, alpha) {
  sumX <- sum(myData)
  prodX <- prod(myData)
  n <- length(myData)
  estimate <- (1/((gamma(alpha^n))*((sumX/(n*alpha))^(n*alpha))))*((prodX)^(alpha-1))*(exp(1)^(-n*alpha)) 
  return(-1*estimate)
} 

Now to optimize, I attempted three different functions from R:

optim(par=0, fn=objFunction, method = "Brent", lower = 0, upper = 10, alpha=2)
nlm(objFunction, momAlpha, myData=myData)
optimize(f=objFunction, c(0,10), alpha=2, maximum=TRUE)

The variable momAlpha, is basically the method of moments estimator for the Alpha, as that would be a good start. Just for completeness:

momAlpha <- (mean(myData)^2)/var(myData)
momBeta <- var(myData)/mean(myData)

These are available in many online references.

Now when I ran the optimization functions above, my results were not clear to me and I need some help understanding:

optim(par=0, fn=objFunction, method = "Brent", lower = 0, upper = 10, alpha=2)
$par
[1] 0.000000005349424
$value
[1] -101196146
$counts
function gradient 
NA       NA
$convergence
[1] 0
$message
NULL

Why is this estimate way out of range?

nlm(objFunction, momAlpha, myData=myData)
$minimum
[1] 0
$estimate
[1] 1.919078
$gradient
[1] 0
$code
[1] 1
$iterations
[1] 0
Warning messages:
1: In f(x, ...) : value out of range in 'gammafn'
2: In f(x, ...) : value out of range in 'gammafn'
3: In f(x, ...) : value out of range in 'gammafn'

The estimate here is nothing but the starting point I provided, why?

optimize(f=objFunction, c(0,2), alpha=2, maximum=TRUE)
$maximum
[1] 1.999934
$objective
[1] -0.2706795

Is this even right?

I am still developing my intuition for the subject, but it seems that I am either doing something wrong, i.e. my objective function is incorrect or the parameter settings of the functions is incorrect or I simply don't understand the way the functions work. I appreciate any help in guiding through this!

Answer 1

You can compute MLE for the gamma distribution using the dglm package, which is available from the CRAN repository. Here is an example run. Note that the two parameters being estimated in this example are the log-mean, which is $\log(\alpha\beta)$, and the log-dispersion, which is $-\log(\alpha)$.

> shape <- 2
> scale <- 1.5
> set.seed(123456)
> myData <- rgamma(n=1000, shape=shape, scale=scale)
> library(dglm)
> fit <- dglm(myData~1, family=Gamma(link="log"), mustart=mean(myData))
> summary(fit)

Call: dglm(formula = myData ~ 1, family = Gamma(link = "log"), mustart = mean(myData))

Mean Coefficients:
            Estimate Std. Error  t value      Pr(>|t|)
(Intercept) 1.117289 0.02197604 50.84124 2.286536e-279
(Dispersion Parameters for Gamma family estimated as below )

    Scaled Null Deviance: 1080.046 on 999 degrees of freedom
Scaled Residual Deviance: 1080.046 on 999 degrees of freedom

Dispersion Coefficients:
              Estimate Std. Error   z value     Pr(>|z|)
(Intercept) -0.7113062 0.04157827 -17.10764 1.301602e-65
(Dispersion parameter for Digamma family taken to be 2 )

    Scaled Null Deviance: 1323.43 on 999 degrees of freedom
Scaled Residual Deviance: 1326.495 on 999 degrees of freedom

Minus Twice the Log-Likelihood: 3992.104 
Number of Alternating Iterations: 2 
> mu <- exp(1.117289)
> shape <- exp(0.7113062)
> scale <- mu/shape
> c(shape, scale)
[1] 2.036650 1.500777

The last line of output gives the MLE for the shape $\alpha$ and the scale $\beta$.

The dglm function is intended to fit mean-dispersion models with link-linear predictors for both the mean and the dispersion of a generalized linear model. The two parameter gamma distribution is a simple special case.

The function uses separate Fisher scoring algorithms for the mean and dispersion parameters, alternating between one iteration of each. For this data, the algorithms converged in two iterations.

Answer 2

The following function, mlgamma() is based on the answer by Gordon Smyth. The function returns the estimated shape and scale of the Gamma function given a vector of sample values. It requires the dglm package to have been installed (install.packages("dglm")) and loaded (library(dglm)).

mlgamma <- function(x) {
  fit <- dglm(
    x~1, 
    family=Gamma(link="log"), 
    mustart=mean(x)
  )
  mu <- exp(fit$coefficients)
  shape <- exp(-fit$dispersion.fit$coefficients)
  scale <- mu/shape
  result <- c(shape, scale)
  names(result) <- c("shape", "scale")
  result
}

An example of its use:

> set.seed(42)
> x <- rgamma(1000, shape=3.6, scale=1.2)
> mlgamma(x)
   shape    scale 
3.548779 1.230362