# Parameter estimation of Gamma Distribution using R

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  0.000000005349424$value
 -101196146
$counts function gradient NA NA$convergence
 0
$message NULL  Why is this estimate way out of range? nlm(objFunction, momAlpha, myData=myData)$minimum
 0
$estimate  1.919078$gradient
 0
$code  1$iterations
 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.999934$objective
 -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!

• Are you asking how to do the ML estimation for the gamma distribution, or are you primarily trying to understand the optimization algorithms that are provided in R? Is the gamma distribution your primary concern or is just a randomly chosen example? The title of your question is too broad to have a succinct answer. Feb 19, 2017 at 23:29
• A bit of both tbh, I am trying to do the ML estimation for the alpha parameter in the gamma distribution, but want to know how to do it in R using the optimization functions in R Feb 19, 2017 at 23:34
• Optimize the log likelihood, not the likelihood itself. You are getting outlandish results for a host of reasons, but by using the log you would have some chance of getting reasonable ones even if you don't correct the other problems.
– whuber
Feb 19, 2017 at 23:35
• @whuber Am I using the optimization functions correctly? Feb 19, 2017 at 23:37
• Not if you're multiplying lots of values together! Trying to optimize such an objective function, even with double-precision floats, is precarious at best: during the search it's almost certain to underflow. Use the log likelihood. Test your objective function to make sure it is correct. Only then try optimization. With experience, you will learn to reparameterize the distribution family so that the parameters are nearly orthogonal, rather than highly correlated, and so that you can optimize either without constraints or with very simple constraints.
– whuber
Feb 19, 2017 at 23:41

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)
 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.

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,
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