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!
Best Answer
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)$.
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.